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PREFACE 


The  study  reported  herein  was  authorized  as  a  part  of  the  Civil  Works  Re¬ 
search  and  Development  Program  by  the  Office,  Chief  of  Engineers  (OCE),  US  Army 
This  particular  work  unit,  Erosion  Control  of  Scour  During  Construction,  is 
part  of  the  Improvement  of  Operations  and  Maintenance  Techniques  (IOMT)  Pro¬ 
gram.  Mr.  James  L.  Gottesman  was  the  OCE  Technical  Monitor  for  the  IOMT  Pro¬ 
gram  during  preparation  and  publication  of  this  report. 

This  study  was  conducted  during  the  period  1  January  1981  through 
31  March  1982  by  personnel  of  the  Hydraulics  Laboratory  of  the  US  Army  Engi¬ 
neer  Waterways  Experiment  Station  (WES)  under  the  general  supervision  of 
Messrs.  H.  B.  Simmons,  Chief  of  the  Hydraulics  Laboratory;  F.  A.  Herrmann,  Jr., 
Assistant  Chief  of  the  Hydraulics  Laboratory;  R.  A.  Sager,  Chief  of  the  Estu¬ 
aries  Division  and  IOMT  Program  Manager;  Dr.  R.  W.  Whalin  and  Mr.  C.  E.  Chatham 
former  and  acting  Chiefs  of  the  Wave  Dynamics  Division,  respectively;  Mr.  D.  D. 
Davidson,  Chief  of  the  Wave  Research  Branch;  and  Dr.  J.  R.  Houston,  Research 
Hydraulic  Engineer  and  Principal  Investigator  for  the  Erosion  Control  of  Scour 
During  Construction  work  unit.  The  Wave  Dynamics  Division  was  transferred  to 
the  Coastal  Engineering  Research  Center  (CERC)  of  the  WES  on  1  July  1983  under 
the  direction  of  Dr.  Whalin,  Chief  of  the  Coastal  Engineering  Research  Center. 
Computer  programming  for  this  study  was  performed  by  Ms.  L.  Chou,  Mathemati¬ 
cian.  This  report  was  prepared  by  Dr.  Houston  and  Ms.  Chou. 

Commanders  and  Directors  of  WES  during  the  conduct  of  this  investigation 
and  the  preparation  and  publication  of  this  report  were  COL  Nelson  P.  Conover, 
CE,  and  COL  Tilford  C.  Creel,  CE .  Technical  Director  was  Mr.  F.  R.  Brown. 
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US  customary  units  of  measurement  used  in  this  report  can  be  converted  to 
metric  (SI)  units  as  follows: 


Multiply _  _ By _  _ To  Obtain 


feet 

0.3048 

metres 

feet  per  second 

0.3048 

metres  per  second 

feet-feet  per  second 

0.0929 

metres-metres  per  second 

feet  per  second  per  second 

0.3048 

metres  per  second  per  second 

inches 

25.5 

millimetres 

square  feet 

0.09290304 

square  metres 
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EROSION  CONTROL  OF  SCOUR  DURING  CONSTRUCTION 


FINITE--A  NUMERICAL  MODEL  FOR  COMBINED 
REFRACTION  AND  DIFFRACTION  OF  WAVES 


PART  I:  INTRODUCTION 


Statement  of  the  Problem 


1 .  It  is  f requen 
improvement  in  the  surf 
tional  beaches,  and  nav 
quarried  rock  or  precas 
When  these  major  struct 
isting  currents.  Shall 
structure  will  cause  bo 
region  by  longshore  or 
is  often  not  compensate 
is  a  scour  hole,  or  ero 
structure.  — 

2.  In  order  to  ensure  structural  stability  and  functional  adequacy  of 
the  works  of  improvement,  any  scour  area  must  be  filled  with  nonerodible  mate¬ 
rial  (sufficiently  stable  to  withstand  the  environmental  forces  to  which  it 
will  be  subjected).  This  may  result  in  additional  quantities  of  material 
being  required  during  construction  that  can  potentially  lead  to  substantial 
cost  overruns.  To  minimize  potential  cost  increases  due  to  scour  during  con¬ 
struction,  it  is  necessary  to  quantify  the  probability  and  ultimate  extent  of 
potential  scour  during  the  scheduled  construction  period.  This  is  an  ex¬ 
tremely  complex  problem  and  quantification  of  the  probability  of  potential 
scour  will  likely  be  site-specific. 

3.  Since  waves  and  wave-induced  currents  play  an  important  role  in  pro¬ 
ducing  scour  near  structures,  it  is  important  to  predict  the  wave  field  in  the 
vicinity  of  structures.  Structures  of  interest  may  have  complex  shapes  and 
the  surrounding  bathymetry  may  be  very  complex.  Both  long-  and  short-period 
waves  mav  attack  the  structure  from  martv  different  angles. 


tly  necessary  to  construct  large  engineering  works  of 


and  nearshore  zon^  to  protect  harbor  entrances,  r^carea- 
igation  channels .  f>  These  structures,  usually  built  from 


t  concrete,  are  placed  in  position  by  crane  or  barge, 
ures  are  erected  in  the  coastal  zone  they  alter  the  ex- 
ow-water  surface  gravity  waves  breaking  on  the  new 


ttom  material  to  be  suspended  and  transported  from  the 


other  currents  that  may  exist  This  removal  of  material 


d  by  an  influx  of  additional  material V—»»d^the  result 
sion,  that  usually  develops  along  the  toe  of  the 


Purpose  of  the  Study  0  1 

'  ■  ■'  ^ 

4.  The  purpose  of  this  study  is  to  develop  a  numerical  model  that  can 
determine  wave  fields  near  arbitrary  structures  in  a  region  of  arbitrary  ami 
variable  bathymetry.  The  model  must  be  capable  of  simulating  both  long  and 
short  waves  approaching  structures  from  any  arbitrary  direction.  In  addition 
the  model  must  be  sufficiently  efficient  to  allow  application  to  real  coastal 


engineering  problems. 


PART  II:  NUMERICAL  MODEL  (FINITE) 


Background 


5.  Short  waves  have  wavelengths  that  are  sufficiently  short  (wavelength 
to  depth  ratio  less  than  approximately  20)  that  propagation  speeds  are  a  func¬ 
tion  of  wave  frequency.  In  the  nearshore  region,  short  waves  with  periods 
from  a  few  seconds  to  approximately  20  sec  play  an  important  role  in  the  move¬ 
ment  of  sediment  and  the  stability  of  coastal  structures. 

6.  Equations  that  govern  long-wave  propagation  over  variable  depths  have 
been  known  for  some  time  (Lamb  1932).  However,  for  shorter  period  waves  where 
frequency  dispersion  is  important,  the  theory  has  included  only  constant 
depths.  Attempts  have  been  made  (Pierson  1951;  Eckart  1952)  to  develop  a  two- 
dimensional  equation  for  combined  refraction  and  diffraction  that  would  govern 
short-wave  propagation  in  a  region  of  variable  depth;  however,  the  equations 
developed  do  not  reduce  to  the  appropriate  simple  refraction  equation  after 
neglecting  the  curvature  of  the  amplitude  function  and  they  also  do  not  reduce 
to  the  linear  long-wave  equation  in  the  case  of  small  water  depth.  In  recent 
years,  Boussinesq-type  equations  also  have  been  used  to  study  the  propagation 
of  short-period  waves.  These  equations  include  terms  that  govern  both  fre¬ 
quency  dispersion  and  nonlinearity.  Two-dimensional  modeling  of  these  equa¬ 
tions  is  difficult  since  the  frequency  dispersion  term  is  third  order,  there¬ 
fore  requiring  a  third-order  accurate  numerical  scheme.  Abbott,  Petersen,  and 
Skovgaard  (1978)  have  presented  a  time-marching,  two-dimensional,  finite  dif¬ 
ference  numerical  model  with  a  third-order  accurate  implicit  solution  scheme. 
However,  a  time-marching  scheme  requires  significant  computational  time.  In 
addition,  the  finite  difference  method  does  not  realistically  model  a  complex 
land-water  interface  as  a  result  of  the  stair-step  representation  and  has  dif¬ 
ficulties  in  allowing  scattered  waves  to  propagate  out  of  a  finite  extent  grid, 
since  these  scattered  waves  are  not  known  a  priori. 

7.  Berkhoff  (1972)  and  Schonfeld  (1972)  have  derived  a  two-dimensional 
wave  equation  that  governs  short-wave  propagation  over  moderately  varying 
depths.  Smith  and  Sprinks  (1975)  gave  a  formal  derivation  of  this  equation. 

In  the  present  study,  a  hybrid  finite  element  model  (FINITE)  is  used  to  solve 
this  wave  equation.  The  method  combines  a  finite  element  solution  over  a 
finite  extent  region  of  variable  depth  with  an  analytical  solution  for  a 


surrounding  infinite  region  of  constant  depth. 


Wave  Equation 


9.  The  propagation  of  periodic,  small  amplitude,  surface  gravity  waves 
over  a  variable  depth  seabed  of  mild  slope  is  governed  by  the  following 
equation . 


V 


(cc^V4>)  +  uj  $  =  0 


(1) 


where 

V  =  horizontal  gradient  operator,  dimensionless 

1/2 

c  =  phase  velocity,  ft/sec*  =  (g/k  tanh  kh) 

2 

g  =  gravitational  constant,  32.174  ft/sec 
k  =  wave  number,  2n/L,  1/ft 
h  =  still-water  depth,  ft 
c^  =  group  velocity,  ft/sec  =  1/2  c(l  +  G) 

G  =  2k  h/sinh  2kh,  ft/sec 

2 

<t>  =  velocity  potential  defined  by  u  =  V<f>  ,  ft  /sec 
u  =  two-dimensional  velocity  vector,  ft/sec 
u)  =  angular  frequency,  2n/T  ,  1/sec 

10.  Equation  1  was  derived  by  Berkhoff  (1972)  and  Schonfeld  (1972)  and 
is  discussed  in  detail  by  Jonsson  and  Brink-Kjaer  (1973).  This  equation  gov¬ 
erns  both  refraction  and  diffraction.  It  reduces  to  the  well-known  "eikonal" 
equation  governing  refraction  by  neglect  of  the  variation  of  the  amplitude 
function  in  the  horizontal  plane.  The  equation  reduces  to  the  diffraction 
Helmholtz  equation  in  deep  or  constant-depth  water  and  to  the  linear  long-wave 
equation  in  shallow  water. 

11.  Finite  element  numerical  models  have  been  used  to  solve  Equation  1. 
Berkhoff  (1972,  1976)  linked  a  finite  element  solution  of  Equation  1  over  a 
variable  depth  region  to  a  source  distribution  for  a  constant-depth  outer 
region.  However,  as  noted  by  Chen  and  Mei  (1974),  Berkhoff  (1972)  did  not  use 
a  proper  functional  with  the  consequence  that  his  global  stiffness  matrix  was 
nonsymmetric  and  thus  inconvenient  numerically  for  all  but  the  simplest  prob¬ 
lems.  Bettess  and  Zienkiewicz  (1977)  also  developed  a  finite  element  solution 


*  A  table  of  factors  for  converting  US  customary  units  of  measurements  to 
metric  (SI)  units  is  presented  on  page  3. 


of  Equation  1.  However,  they  used  infinite  elements  to  cover  the  constant- 
depth  outer  region.  The  shape  function  used  in  the  infinite  elements  had  an 
exponentially  decreasing  term  in  the  direction  away  from  the  inner  region. 

The  choice  of  a  decay  length  was  somewhat  arbitrary;  but  the  solutions  were 
not,  in  general,  too  sensitive  to  the  exact  value.  A  disadvantage  of  this 
technique  is  that  the  infinite  elements  increase  the  number  of  equations  to  be 
solved.  If  a  problem  requires  a  large  number  of  elements  because  the  region 
of  interest  is  large  and  the  incident  wavelengths  are  short,  this  solution  can 
require  substantial  computational  time. 

12.  Equation  1  also  has  been  solved  by  a  parabolic  approximation  (Radder 
1979).  The  approximation  is  derived  from  splitting  the  wave  field  into  trans¬ 
mitted  and  reflected  components  and  then  neglecting  the  reflected  components. 
This  approach  is  applicable  to  some  propagation  problems  but  is  not  appropri¬ 
ate  for  problems  involving  wave  interaction  with  coastal  structures  such  as 


breakwaters . 


Finite  Element  Solution 


13.  Equation  1  is  solved  by  program  FINITE  using  a  hybrid  finite  element 
method  originally  developed  by  Chen  and  Mei  (1974)  to  solve  the  diffraction 
Helmholtz  equation  in  a  constant-depth  region.  Space  is  divided  into  two  re¬ 
gions  as  shown  in  Figure  1  (finite  inner  region  A  and  infinite  outer  re¬ 
gion  B) .  Conventional  finite  elements  are  used  in  the  variable  depth  re¬ 
gion  A  .  A  single  superelement  is  used  to  cover  the  constant-depth  infinite 
region  B  .  Variational  principles  are  used  that  incorporate  the  matching  con¬ 
ditions  between  the  regular  elements  and  the  superelement  as  natural  condi¬ 
tions.  Thus  a  symmetric  global  stiffness  matrix  is  obtained  that  is  very  ad¬ 
vantageous  for  highly  complex  problems. 

14.  The  variational  principle  for  the  boundary  value  problem  requires 
that  the  following  functional  be  stationary  with  respect  to  arbitrary  first 
variation  of  the  velocity  potential  4>  : 


2 

<p  dA 


J  ^  c  -  ^ 


cc  ad). 


■  *  V  '  .< 


"  *  *  M  *  *  **«"»'* 


-  boundary  b 

(AT  INFINITY) 


where  4>^  and 


REGION  S 


Figure  1.  Regions  of  computation 


are  the  incident  and  reflected  wave  velocity  potentials, 


respectively;  n  is  a  unit  normal;  and  the  last  two  integrals  are  line  inte 
grals  at  infinity.  Analogous  to  the  derivation  of  Chen  and  Mei  (1974),  this 
functional  can  be  rewritten  as  follows: 


-fflt 
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ac<t>B  -  ♦j) 


frl*m 


Sin] 


where  <J>g  and  <J>^  are  the  velocity  potentials  in  regions  B  and  A  ,  re¬ 
spectively,  and  n^  is  a  unit  normal  to  the  boundary  separating  regions  B 
and  A  . 

15.  Note  that  all  integrals  are  evaluated  within  region  A  or  along 
3A  .  Thus  the  variational  principle  is  a  ’realized  one.  As  discussed  by 
Aranha,  Mei,  and  Yue  (1979),  this  variational  principle  can  be  replaced  by  an 
equivalent  weak  formulation. 

16.  The  inner  region  A  is  assumed  to  have  a  variable  depth  and  to  be 
of  finite  extent.  This  region  is  subdivided  into  finite  elements.  Here  the 
elements  are  triangular  with  simple  linear  shape  functions.  The  infinite 
region  B  (or  semi-infinite  if  a  harbor  located  along  an  infinite  coastline 
is  of  concern)  is  assumed  to  have  a  constant  depth  and  is  covered  with  a 
single  superelement.  Since  region  B  has  a  constant  depth,  the  governing 
equation  is  the  diffraction  Helmholtz  equation.  An  analytical  solution  for 
the  velocity  potential  in  region  B  is  well  known  and  can  be  expressed  as 
follows : 

00 

<t>D  =  H  (kr)(a  cos  n©  +  p  sin  n0)  (4) 

B  n  n  Kn 

n=0 

where  a  and  p  are  constant  and  unknown  coefficients,  H  (kr)  is  Hankel 
n  rn  n 

function  of  the  first  kind,  and  r  and  0  are  radial  and  angular  variables 
in  polar  coordinates.  For  a  semi-infinite  region  B  and  a  straight  infinite 
coastline,  <J>D  can  be  expressed  as  follows: 

D 

00 

<J>D  =  Y'  a  H  (kr)  cos  n0  (5) 

B  n  n 

n=0 

The  velocity  potentials  given  in  Equations  4  and  5  satisfy  the  Sommerfeld 
radiation  condition  that  the  scattered  waves  must  behave  as  outgoing  waves  at 
infinity.  Thus  region  B  can  be  considered  to  be  a  single  superelement  with 
shape  function  given  by  Equations  4  and  5. 

17.  If  the  shape  functions  are  used  to  evaluate  the  integrals  of  Equa¬ 
tion  3  and  the  functional  is  extremized  with  respect  to  the  unknowns,  a  set  of 
linear  algebraic  equations  is  obtained.  Of  course,  the  infinite  series  given 
by  Equations  4  and  5  must  be  truncated  at  some  finite  extent.  The  number  of 
terms  that  must  be  retained  depends  upon  the  incident  wavelength  and  may  be 
found  by  increasing  the  number  of  terms  until  the  solution  is  insensitive  to 


the  addition  of  further  terms.  Solution  of  the  boundary  value  problem  thus 
reduces  to  the  solution  of  N  linear  algebraic  equations  for  N  unknowns 
(where  N  is  the  number  of  node  points  in  the  finite  element  discretization 
plus  the  number  of  unknowns  in  the  truncated  series).  That  is, 

[K]  {<!<}  =  {Q}  (6) 

N  x  N  N  x  1  N  x  1 

The  symmetric  complex  coefficient  matrix  [K]  is  in  general  large,  sparse, 
and  banded.  It  can  be  stored  and  manipulated  in  the  computer  in  a  packed  form 
(Chen  and  Mei  1974).  The  packed  form  is  chosen  to  be  a  rectangular  array  (N 
variables  in  length  and  the  semibandwidth  in  width).  Only  elements  of  [K] 
on  and  above  the  diagonal  and  within  the  band  width  need  to  be  stored  in  the 
packed  form. 

18.  Although  the  packed  form  of  [K]  greatly  reduces  the  required  com¬ 
puter  memory,  the  problems  discussed  later  are  so  large  that  even  memory  re¬ 
quirements  of  the  packed  form  of  [K]  are  excessive.  However,  since  the  sym¬ 
metric  coefficient  matrix  is  positive  definite,  a  solution  is  possible  by 
elimination  methods  without  pivoting.  Without  pivoting,  elimination  performed 
using  one  row  affects  only  the  triangle  of  elements  within  the  band  below  that 
row.  Thus  the  packed  form  of  [K]  can  be  partitioned  into  several  smaller 
blocks.  Using  Gaussian  elimination,  only  two  blocks  at  a  time  are  involved  in 
the  reduction  and  back  substitution  with  the  remainder  of  the  blocks  kept  in 
peripheral  storage.  This  technique  allows  the  solution  of  extremely  large 
matrices . 


Verification 

19.  To  verify  this  numerical  model,  comparisons  were  made  between  the 
finite  element  calculations  and  both  an  analytical  and  a  numerical  solution 
for  the  interaction  of  waves  with  a  circular  island  on  a  paraboloidal  shoal. 
Figure  2  is  a  sketch  of  the  problem.  Hom-ma  (1950)  presented  the  analytical 
solution  to  the  long-wave  equation  for  plane  waves  incident  upon  this  island. 

20.  Many  investigators,  including  Berkhoff  (1972,  1976)  and  Bettess 
and  Zienkiewicz  (1977),  have  compared  their  finite  element  solutions  of  the 
long-wave  equation  with  Hom-ma ' s  solution.  For  example,  Berkhoff  (1972,  1976) 
used  156  elements  to  cover  the  paraboloidal  shoal  (actually  one-half  of  the 
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Figure  2.  Circular  island  on  a  paraboloidal  shoal 


shoal  since  by  symmetry  the  solution  required  calculations  for  just  one-half 
of  the  shoal).  This  number  of  elements  was  not  sufficient  to  adequately  re¬ 
solve  features  of  the  problem;  consequently,  Berkhoff's  finite  element  solu¬ 
tion  differed  from  Hom-ma ' s  analytical  solution  by  as  much  as  70  percent  at 
one  location.  Since  Berkhoff  did  not  attempt  to  use  a  finer  grid,  it  is 
likely  that  his  coarse  grid  required  a  moderate  computational  time.  The 
model  described  in  this  paper  would  require  less  than  0.1  sec  of  computa¬ 
tional  time  on  a  CRAY-1  computer  to  solve  Equation  6  using  Berkhoff's  grid. 
Bettess  and  Zienkiewicz  (1977)  also  used  a  fairly  coarse  grid  similar  to 
Berkhoff’s  grid. 

21.  Lautenbacher  (1970)  used  an  integral  equation  solution  to  solve  the 
long-wave  equation  for  waves  interacting  with  a  circular  island  on  a  shoal 
with  linear  side  slopes.  He  used  a  coarse  circular  mesh  grid  with  only 

130  points.  Because  the  resulting  coefficient  matrix  was  full,  the  computa¬ 
tional  time  required  for  a  solution  was  60  min  on  an  IBM  7094  computer. 

22.  Figure  3  shows  a  finite  element  grid  with  2,640  elements  used  by 
the  model  described  in  this  report  to  solve  the  problem  of  the  interaction  of 
long  waves  with  a  circular  island  on  a  paraboloidal  shoal  (by  symmetry  only 


Figure  3.  Finite  element  grid  for  paraboloidal  island 

(2,640  elements) 


half  the  shoal  needs  to  be  considered).  Figure  4  shows  comparisons  between 
Hom-ma ' s  analytical  solution  and  the  finite  element  model  solution  for  incident 
waves  with  five  different  periods.  The  agreement  is  excellent  with  only  slight 
differences  for  the  240-sec  waves  (resulting  from  lower  resolution  of  the  inci¬ 
dent  wave  for  shorter  period  waves).  The  computational  time  for  solution  is 
approximately  4  sec  on  a  CRAY-1  computer. 
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Figure  4.  Solutions  without  dispersion  (2,640  element  grid) 


23.  Jonsson,  Skovgaard,  and  Brink-Kjaer  (1976)  show  that  for  a  wave 
with  a  240-sec  period  interacting  with  the  circular  island  on  a  paraboloidal 
shoal  the  effect  of  frequency  dispersion  is  not  particularly  significant.  The 
ratio  of  wavelength  to  water  depth  for  this  case  is  approximately  11.  However 
for  a  120-sec  incident  wave  (wavelength  to  water  depth  ratio  of  less  than  5), 
frequency  dispersion  is  quite  significant.  In  order  to  maintain  a  resolution 
of  a  120-sec  wave  that  is  approximately  equal  to  that  obtained  for  the  240-sec 
wave  using  the  2,640-element  grid,  it  is  necessary  to  reduce  element  side 
lengths  by  a  factor  of  approximately  2.  This  reduction  results  in  a  quadru¬ 
pling  of  the  number  of  elements.  Figure  5  shows  a  finite  element  grid  with 
10,560  elements  used  to  calculate  the  interaction  of  a  120-sec  wave  with  the 
island . 


Figure  5.  Finite  element  grid  for  paraboloidal  island 

(10,560  elements) 

24.  Figure  6  shows  a  comparison  between  the  analytical  solution  of  the 
long-wave  equation  by  Hom-ma  (1950)  and  the  finite  element  solution  (grid  of 
Figure  5)  for  a  120-sec  incident  wave.  Figure  7  also  shows  a  comparison  be¬ 
tween  a  numerical  solution  of  Equation  1  using  an  orthogonal  collocation  solu¬ 
tion  and  the  finite  element  model  solution  of  this  report.  In  both  cases 
there  is  excellent  agreement.  The  effect  of  frequency  dispersion  is  quite  sig¬ 
nificant  as  shown  in  Figure  8,  where  the  solutions  of  the  long-wave  equation 
and  Equation  1  are  overlapped.  In  fact,  the  inclusion  of  frequency  dispersion 
is  much  more  significant  a  factor  than  is  resolution  of  the  wave  form  for  this 
particular  case.  For  example,  Figure  9  shows  a  solution  including  dispersion 
using  the  2,640-element  grid.  The  agreement  between  the  finite  element  and 
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Figure  6.  Solutions  without  dispersion  (10 ,560-element  grid) 

-  ORTHOGONAL  COLLOCATION  SOLUTION  (JONSSON  ETAL.  1976) 
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Figure  8.  Effect  of  dispersion 

colloc.it  ion  solutions  is  much  better  using  the  2,640-element  grid  and  including 
dispersion  than  is  obtained  using  the  10,560-element  grid  without  dispersion. 

25.  The  computational  time  requirement  to  solve  a  problem  using  the 
1 0 , 560- e 1 ement  grid  is  less  than  1  min  on  a  CRAY-1.  This  is  extremely  modest 
considering  the  very  large  size  of  this  finite  element  grid.  The  computa¬ 
tion.!  i  time  of  the  model  is  proportional  to  the  number  of  nodes  times  the 
bandwidth  of  the  coefficient  matrix  squared.  This  grid  has  more  than  50  times 
the  number  of  nodes  contained  in  the  grid  used  by  Berkhoff  (1972,  1976)  and  a 
bandwidth  approximately  3.3  times  greater.  Thus  the  computational  time  re¬ 
quirement  :s  almost  600  times  greater  for  this  grid  than  for  Berkhoff 's  grid. 
The  computer  memory  requirements  also  are  very  modest  as  a  result  of  the  parti¬ 
tioning  of  the  coefficient  matrix.  Although  the  packed  form  of  the  coefficient 
matrix  had  approximately  1.7  million  terms  for  this  problem,  only  45,000  terms 
we i e  in  central  memory  at  any  given  time  with  the  remainder  in  peripheral 


storage . 
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Figure  9.  Solutions  with  dispersion  (2,640-element  grid) 


Comparison  with  Three-Dimensional  Numerical  Model 


26.  Equation  1  was  derived  assuming  that  the  bathymetry  was  slowly  vary 
ing.  This  mild  slope  approximation  is  valid  for  the  slow  bathymetric  varia¬ 
tions  of  the  paraboloidal  shoal  shown  in  Figure  2.  However,  the  approximation 
is  not  appropriate  for  many  practical  problems  involving  the  interaction  of 
waves  with  man-made  structures.  For  these  problems,  it  is  necessary  to  use  a 
fully  three-dimensional  model  such  as  the  hybrid  three-dimensional  finite  ele¬ 
ment  developed  by  Yue,  Chen,  and  Mei  (1976).  Of  course,  a  three-dimensional 
model  requires  substantial  computational  time.  Therefore  it  is  of  interest  to 
consider  how  well  a  two-dimensional  mild-slope  model  compares  with  a  three- 
dimensional  model  for  a  problem  where  the  mild-slope  approximation  is  strongly 
violated . 

27.  Yue,  Chen,  and  Mei  (1976)  used  a  three-dimensional  model  to  con¬ 
sider  the  interaction  of  small  amplitude  plane  waves  with  an  elliptic  island 
on  a  circular  base.  Figure  10  illustrates  thp  problem  and  shows  the  finite 


Figure  10.  Elliptic  island  and  three-dimensional  grid 

element  grid  used  by  Yue,  Chen,  and  Mei  (1976).  Bottom  side  slopes  vary  from 

IV  on  1H  to  IV  on  2H  and  violate  the  mild-slope  assumption. 

28.  Figure  11  shows  the  two-dimensional  finite  element  grid  used  by  the 

model  described  in  this  report  to  simulate  the  interaction  of  small  amplitude 

plane  waves  with  the  elliptic  island.  Figures  12  and  13  show  the  interaction 

of  waves  (k  a  =  1  and  two  incident  directions)  with  the  elliptical  island, 
o 

The  amplification  factors  and  phases  around  the  elliptical  island  calculated 
by  the  three-dimensional  model  of  Yue,  Chen,  and  Mei  (1976)  and  the 


Figure  11.  Finite  element  grid  for  elliptical  island 

two-dimensional  model  described  in  this  report  are  shown.  The  maximum  differ¬ 
ence  in  amplification  factors  is  not  much  greater  than  10  percent.  Similar 
agreement  was  found  for  greater  values  of  k  a 

29.  Although  the  two-dimensional  model  cannot  perfectly  reproduce  the 
three-dimensional  model  results,  the  difference  may  be  within  the  accuracy 
requirements  of  many  engineering  applications.  Of  course,  the  computational 
requirements  of  the  two-dimensional  model  are  very  modest  compared  with  the  re 
quirements  of  the  three-dimensional  model.  For  the  problem  of  waves  interact¬ 
ing  with  the  elliptical  island,  the  three-dimensional  model  required  2.6  min 
of  computational  time  on  an  IBM  370/168  computer  for  each  wave  period.  The 


two-dimensional  model  required  less  than  2  sec  of  computational  time  on  a 
CRAY-1  computer.  Part  of  the  difference  in  computational  time  is  a  conse¬ 
quence  of  the  unknown  relative  speeds  of  the  two  computers. 

30.  The  difference  in  computational  time  of  the  two-  and  three- 
dimensional  models  becomes  more  significant  for  larger  problems.  For  example, 
Yue,  Chen,  and  Mei  (1976)  estimate  that  approximately  6  hr  of  computational 
time  on  an  IBM  370/168  computer  would  be  necessary  to  calculate  the  interac¬ 
tion  of  a  wave  with  a  period  of  8.5  sec  or  greater  with  a  particular  offshore 
harbor.  If  a  resolution  of  10  grid  points  per  wavelength  is  maintained,  the 
two-dimensional  model  would  require  an  estimated  computational  time  of  only 
0.5  to  1  min  to  perform  the  same  calculation  on  a  CRAY-1  computer. 


Comparison  with  Laboratory  Experiments 

31.  Putnam  and  Arthur  (1948)  performed  the  pioneering  laboratory  experi¬ 
ments  that  considered  diffracted  waves  in  the  lee  of  an  impermeable  breakwater 
located  in  constant-depth  water.  Mobarek  (1962)  performed  similar  experiments 

except  that  a  linear  bottom  slope  was  used  in  the  lee  of  the  breakwater  and  a 

2 

very  small  model  (72  ft  )  was  used.  In  order  to  study  the  phenomenon  of  com¬ 
bined  refraction  and  diffraction  near  structures,  laboratory  experiments  were 
performed  by  Hales  (1980).  In  these  experiments  an  impermeable  breakwater  was 
located  perpendicular  to  a  straight  coastline.  A  linear  bottom  slope  extended 
from  some  distance  in  front  of  the  breakwater  to  the  shoreline. 

32.  The  laboratory  facility  used  in  these  tests  covered  an  area  of  ap- 

2 

proximately  2,500  ft  .  The  breakwater  was  15  ft  long  and  1  in.  thick.  The 
water  depth  in  the  facility  decreased  from  1.0  ft  to  0.0  ft  over  a  distance  of 
20  ft.  The  sidewalls  that  laterally  bounded  the  facility  were  curved  to  fol¬ 
low  wave  orthogonal s.  A  plunger-type  wave  maker  was  used  to  generate  small 
amplitude  sinusoidal  waves  that  approached  the  breakwater  at  an  angle.  An 
array  of  32  para  1 lei -wi re ,  resistance-type  sensors  was  used  to  measure  the 
waves.  The  gages  were  supported  by  a  large  stand  so  that  only  the  measuring 
wires  of  the  gages  were  in  the  water  (i.e.,  the  gages  did  not  require  indi¬ 
vidual  feet  for  support).  Information  was  recorded  and  analyzed  by  a 
mini  computer . 

33.  One  problem  in  simulating  these  laboratory  tests  numerically  is  that 
the  waves  break  in  the  hydraulic  model  near  the  shoreline  and  thus  dissipate 


their  energy.  There  is  no  mechanism  to  dissipate  energy  in  the  numerical  model 
described  in  this  paper.  However,  dissipation  can  be  simulated  by  allowing 
waves  to  continue  to  propagate  out  of  the  problem  area.  Figure  14  illustrates 
schematically  how  this  is  done.  The  breakwater  and  the  linear  slope  are 
numerically  modeled  only  to  the  point  where  breaking  occurs.  The  depth  is 
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Figure  14.  Wave  interaction  with  impermeable  breakwater 


then  increased  to  the  depth  of  the  semi-infinite  region  surrounding  the  inner 
region  and  the  waves  are  allowed  to  radiate  away  from  the  inner  region.  Fig¬ 
ure  15  shows  the  finite  element  grid  used  for  this  simulation. 

34.  Figure  16  shows  a  typical  comparison  between  the  laboratory  mea¬ 
surements  and  the  finite  element  calculations.  Also  shown  is  a  uniformly 
valid  asymptotic  solution  derived  recently  by  Liu  and  Lozano  (1979).  The 
solution  derived  by  Liu  and  Lozano  (1979)  is  in  excellent  agreement  with  the 
laboratory  tests.  The  finite  element  calculations  agree  quite  well  in  the 
shadow  zone  with  the  laboratory  tests.  The  agreement  is  not  as  good  outside 
the  shadow  zone.  The  difference  probably  is  attributable  to  the  artificial 
increase  in  depth  to  allow  the  waves  to  radiate  from  the  inner  region.  This 
depth  transition  would  cause  some  energy  to  reflect  back  into  the  inner  region 
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Figure  16.  Comparison  with  laboratory  experiments 


PART  III:  INPUT  DATA  PREPARATION  FOR 
FINITE  ELEMENT  SOLUTION 


35.  The  following  example  problem  is  presented  to  demonstrate  the  steps 
and  procedures  for  preparing  a  finite  element  grid  to  be  used  for  the  numerical 
model  solution  of  combined  refraction  and  diffraction.  This  example  outlines 
the  methodology  for  data  manipulation  on  the  computer  facilities  at  the  US  Army 
Engineer  Waterways  Experiment  Station  (WES). 


Data  Preparation  for  Computer  Plotting 


36.  Assume  it  has  been  determined  that  an  area  of  interest  such  as  out¬ 


lined  in  Figure  17  is  required  to  be 


Figure  17.  Example  problem  for 
demonstrating  the  preparation 
of  input  data  for  hybrid 
finite  model 


igitized  for  computer  processing  of  the 
hybrid  finite  element  model.  The 
initial  determinations  after  the  over¬ 
all  region  has  been  delineated  include 
an  orientation  so  that  the  elements  of 
the  grid  can  be  numbered  in  continuous 
fashion.  Also,  the  nodes  of  the  ele¬ 
ments  should  be  numbered  in  a  continu¬ 
ous  fashion  (counterclockwise).  The  x 
and  y  coordinates  of  each  node  of  the 
elements  are  digitized  in  continuous 
fashion  and  recorded  for  verification 
later.  If  the  program  CONVER  is  used 
for  converting  the  measurement  of  digi¬ 
tized  x  and  y  coordinates  to  proto¬ 
type  units,  it  is  necessary  that  all  x 
and  y  values  be  positive.  When  the 
grid  has  been  digitized  on  magnetic 
tape,  the  tape  is  assigned  a  number  for 
relocation  later. 

37.  To  convert  the  digitized 
ASCII  tape  to  BCD  tape,  run  program 
DIG1  (Figure  18).  The  steps  for  per¬ 
forming  such  an  operation  include: 
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a.  Sign  on  DPS1  Honeyvell  System  at  WES.  The  telephone  numbers 
for  both  300  Baud  and  1200  Baud  are  2162,  2171,  3561,  and  3571. 

b.  Type  in 

OLD  ROHH314/DIG1 ,R 

CLEAR 

SAVE  DIG1 

c.  Type  in 
EDIT 

The  system  will  respond  with  a 

d.  Type  in  "PS : /IDENT/"  after  the 

The  system  will  search  and  print  the  line  containing  the  char¬ 
acter  string  "ident"  as: 

20$ : IDENT : USERID .USERNAME 

e.  Type  in 

RVS: /USERID, USERNAME/: /YOUR  USER  ID,  YOUR  NAME/ 

The  system  will  replace  the  old  user's  ID  and  user's  NAME  by 
the  new  user's  ID  and  NAME,  as  desired. 

f.  Type  in 
PS : / 720/ 

The  system  will  respond  with 

720$TAPE9:01,X1D, ,TAPEN0, .ASCII  TAPE  NAME , , DEN8 

g.  Type  in 

RVS ./TAPENO, .ASCII  TAPE  NAME , ,DEN8/ : /ACTUAL  TAPE  NO, .ACTUAL 
ASCII  TAPE  NAME, .DENSITY  OF  THE  TAPE/ 

where  TAPENO  has  been  previously  assigned  to  the  digitized 
magnetic  tape.  ASCII  TAPE  NAME  can  be  (1-8)  arbitrary 
alpha-numeric  characters.  Parameter  DEN8  specified  800 
density  while  DEN16  indicates  1600  density. 

h.  Type  in 
PS:/740/ 

The  system  will  respond  with 

740$ :TAPE9 : 03 ,X3D , , , ,BCD  TAPE  NAME 

i.  Type  in 

RVS: /BCD  TAPE  NAME/ : /ACTUAL  BCD  TAPE  NAME/ 

2-  Type  in 
PS : /750/ 

The  system  will  respond  with 

750$ :MSG2: save  03, USER  NAME .USERID, BCD  TAPE  NAME 

k.  Type  in 

RVS : /USER  NAME .USERID, BCD  TAPE  NAME/ : /YOUR  NAME, YOUR  ID, ACTUAL 
BCD  TAPE  NAME/ 

(Note:  BCD  TAPE  NAME  should  be  the  same  name  as  that  used 
in  Line  740 . ) 

l.  Type  in 
RESAVE  DIG1 

The  system  will  respond  with 
DATA  SAVED  -  DIG1 
End  of  file 
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m.  Type  in 
DONE 

The  system  will  respond  with  a  (CARDIN  MODE)  and  wait  for 
the  next  command. 

n .  Type  in 
RUN 

The  system  will  respond  with 
SNUMB  XXXXXT 

38.  The  BCD  tape  contents  will  be  saved  as  a  permanent  file  by  running 
program  TP2FIL  (Figure  19)  in  the  following  manner. 

a.  Sign  on  the  DPS1  computer  at  WES 

b.  Type  in 

OLD  R0HH314/TP2FIL.R 
CLEAR 

SAVE  TP2FIL 

c .  Type  in 
EDIT 

d.  Type  in 

PS : /IDNET/ 

e.  Replace  USERID  and  USER  NAME  with  actual  user's  ID  and  user's 

NAME 
Type  in 

RVS: /USERID, USER  NAME/ : /ACTUAL  USER’S  ID,  USER'S  NAME/ 

f.  Type  in 
PS:/60/ 

The  system  will  respond  with 

60$ : TAPE9 :F1 ,X1D, ,BCD  TAPENO, ,BCD  TAPE  NAME 

£•  Type  in 

RVS: /BCD  TAPENO,, BCD  TAPE  NAME/ : /ACTUAL  BCD  TAPE  NUMBER, .ACTUAL 
BCD  TAPE  NAME/ 

(Note:  The  BCD  TAPE  NUMBER  is  obtained  from  the  log  file 
of  the  printout  after  program  DIG1  has  been  run,  example 

OLIt  TF2FIL 
♦  LIST 

io»*k. j 

20*  :  IOENTJUSERIDfUSERHAKE 
30*  UTILITY 

40*;FUTIL:F3fF2f  COF'Y/IR/ » HOLP/F2/ 

50*  :FUTIL  :ri  fF2fC0PY/1F/fRW£i/F2/ 

60* : TAPE9 : F 1 f XIDf f BCD  TAPENO »  f  BCDTAF'E  NAME 

7o*:file:fc.:<2s>il 

80*:data:f3i .cofy 

voini ut.asis 

100* : ENDCOPY 

iio*:frogran:tsconv 
120* :file: i*,X2R 

130t:FlLE:0T,  .ILfNEU.PF  NAME 
1 40* ! EHD JOB 


t 

Figure  19.  Program  TP2FIL  to  save 
BCD  tape  as  a  permanent  file 
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shown  in  Figure  20.  The  same  BCD  TAPE  NAME  should  be  used 
as  that  previously  used  in  DIG1.) 

h.  Type  in 
PS:/ 130/ 

The  system  will  respond  with 

130$ :FILE : 0T, 1L,NEW, PF  NAME 

i.  Replace  the  PF  NAME  with  actual  permanent  file  name 
Type  in 

RVS:/PF  NAME/: /ACTUAL  PF  NAME 

(The  same  PF  NAME  as  BCD  TAPE  NAME  is  suggested.) 

j. .  Type  in 

RESAVE  TP2FIL 

k.  Type  in 
DONE 

l.  Type  in 
RUN 

The  system  will  respond  with 
SNUMB  XXXXXT 

OPERATOR  STARTED  WITh  *100*8  TO*  flLC  CODE  OS  CE  600  STL  GC  PHI  100*8  100*8  0001  8*016  OOOHARTIN8C0 

T*UM0?  100*8  WAS  CHAINED  TO  CONTROL  'ASTER  T  A  PE  N A  HE  ROCRHJC 
•  NORMAL  TERMINATION  AT  0 21  S* S  t«5000  S W» 00 00 00 000000 


START 

9.29* 

LINES 

213 

PAOC  0.0018 

i/o  o.ooi 

IU  S 

MEMORY 

STOP 

9.299 

LIMT 

20*80 

LIMIT  0.3000 

LIMIT 

cu  $ 

M  •  T 

SWAP 

0.000 

LAPSE 

0.0U5 

f  c 

0 

TYPE 

SUSY 

1  P/AT  t P/RT 

IS/fC  MS/*! 

ADDRESS 

f  * /PK  * 

3* 

R 

'SU  500 

• 

*2 

*  * 

36  36 

0-12-20 

R  • 

R 

'SUSOO 

• 

1 1 

0  0 

1  1 

0-12-1 3 

OT 

0 

TAPE-9 

T  386 

0/01 

*06  0 

0-16-05 

RU0180 

03 

D 

1  APE-9 

352 

0/03 

*9  0 

0-18-05 

*100*8 

p  • 

STOUT 

L  • 

R 

'SU *50 

• 

1290 

0  0 

1200  1200R 

0-08-02 

• 

•L 

R 

'SU  SCO 

• 

60 

0  0 

S00  50QR 

0-1 2-1 3 

LIST  99  LINES 

A  C -52  11*  LINES 


activity  COST  SUMMARY 

RESOURCE  RESOURCE  BILLING 

DESCRIPTION  USAGE  UNITS 

PROCESSOR  (SECS)  6  10 

DISC  1/3  TIME  (SECS)  1  1 

TAPE  I/O  TIME  (SECS)  2  1 

Figure  20.  Example  showing  how  BCD  TAPE  NUMBER  is  obtained  from 
log  file  of  printout  after  program  DIG1  is  run 

39.  Digitizing  errors  may  be  edited  from  the  tape  with  the  following 
procedures . 

a.  Sign  on  the  DPS1  computer  at  WES 

b.  Type  in 
OLD  PFNAME 

(PFNAME  is  the  permanent  file  name  used  in  TP2FIL.) 


fl  _  If 


to  indicate  that  it  is  in 


c . 


Type  in 
EDIT 


The  system  will  respond  with  a 
the  EDIT  mode. 


d.  Use  the  editor  command  to  correct  all  errors. 

e .  Type  in 
DONE 

The  system  will  respond  with  a  and  return  to  the  CARDIN 
mode  waiting  for  the  next  command. 

f.  Type  in 
RESAVE  PFNAME 

(This  should  be  the  same  file  name  as  PFNAME.) 


40.  Program  CONVER  will  convert  the  digitized  values  of  the  x  and  y 
coordinates  to  prototype  units  (Figure  21). 

a.  Sign  on  the  DPS1  computer  at  WES 

b.  Type  in 

OLD  ROKH314/CONVER ,R 
CLEAR 

SAVE  CONVER 

c.  Type  in 
EDIT 

d.  Replace  USERID  and  USER  NAME  by  actual  user's  ID  and  NAME 
Type  in 

RVS:/USERID, USERNAME/ :/ACTUAL  USER’S  ID,  ACTUAL  USER’S  NAME/ 

e.  Type  in 

PS : /PARAMETER/ 

The  system  will  respond  with  the  statement 

150  PARAMETER  XMO=0 , YMO=0 ,XM1=12 . ,  YM1=0. 

(where  XMO,  YMO,  XM1,  and  YM1  are  coordinates  of  reference 
points  in  the  prototype  units) 

f .  Type  in 

RVS:/XM0=0LD  VALUE , YM0=0LD  VALUE/ : /XM0=NEW  VALUE , YM0=NEW  VALUE/ 
RVS:/XM1=0LD  VALUE , YM 1 =0LD  VALUE/ : /XM0=NEW  VALUE ,YM1=NEW  VALUE/ 

g.  Locate  and  replace  with  parameter  N1  and  N2  in  Line  160 
Type  in 

PS : / 160/ 

RVS:/N1=0LD  VALUE, N2=0LD  VALUE/ : /N1=NEW  VALUE, N2=NEW  VALUE/ 

h.  Locate  and  change  the  format  statement  in  Line  290 
Type  in 

PS : /290/ 

RVS : /FORMAT (4 ( i6 , 2F7 . 2 ) ) / : /THE  DESIRED  FORMAT/ 


OLD  CONUER 
*L  1ST 

100!  IN 

110»: I  DENT (USER ID. USERNAME 
120*:0F  TIONJFORVRAN 
125UUSE  :  .GTLIT 
130*: FORTY IHFORM.NLNO 

150  PARAMETER  X  M  0  -  0 .  ,  V  H  Q  -0.  » X II 1  -  12.  » Y  H 1  -  0  • 

160  PARAMETER  HI  185 1 » N2= 3973 

165  PARAMETER  N3- N2-N141 

170  DINE N SION  XP ( N3 > » YP ( N3 ) « NTT ( H3 ) 

190  DIST-SQRTC <ARS(XH1-XK0) ) **2 .  4  < ARS C YH1- VHO ) )**2.  ) 

192  DO  255  I-I.N3 

19  4  NTT<  D-I  +  Nl-l 

196  255  CONTINUE 

200  REAI .1,100)  XO.YO 

210  READ ( 1 , 100 )  X1,Y1 

220  100  FORMAT  (2X,F4*0»2X,F4.0> 

230  RE AD < 1 , 1 00 >  < XP ( I > , V P < I > , I -1 ,N3 ) 

240  CALL  CONUER ( X0»  YD , X 1 » Y1 «  Nl »H2 . XP , YP, DI ST  »XMD. VHO ) 

270  WRITE <43, 200)  < NTT < I ) , XP < I ) , VP < I > , I = 1 » H3 ) 

275  PRINT  200. < NTT < I ) »XP< I ) »VP< I ) , I =1 iN3> 

290  200  F0RMAT<4< I6.2F7.2) ) 

300  REWIND  1 
310  CALL  DET ACH < 1 , , 1 
320  STOP 
330  END 

340  3URR0UTINE  CONUER < XO , VO , XI , Y 1 , H 1 , N2 , X . Y r DI ST , XKO , VHO ) 
350  DIMENSION  X<1), Y<1> 

360  HYP- 3QRT<  < Y 1 -Y0) **24 < X 1 -XO ) **2 > 

370  SCALE  DIST/HYP 

380  SINE-< <Y1-Y0)*SCAL£)/PIST 

390  COS INE- < (XI -XO) *SCALE)/DIST 

395  N3-N2-N141 

400  DO  10  I-1.N3 

410  XX- < (X<I )-XO)*COSINE+(Y< I ) - YO > *S INE > *SC ALE+XMO 

420  Y  <  I )  -  <-<:<<  I  )-X0)*SINE4< Y<  I  )-YO)*COSINE)*SCALE+YMO 

430  X(I)pXX 

440  10  CONTINUE 

450  RETURN 

460  END 

470*:EXECUTE 

43o*:limits:o2,16k. ,3000 
4V0*:F'RMFl:01  ,R,L,USERID/PF  NAME 
500* 1ENDJ0R 


« 

Figure  21.  Program  CONVER  to  convert  digitized 
values  of  x  and  y  coordinates  to  prototype 

values 


Locate  and  replace  the  parameter  of  USERID  and  PFNAME  in 
Line  490 
Type  in 
PS:/490/ 

RVS: /USERID/: /ACTUAL  USER'S  ID/ 

RVS : /PFNAHE/ : / 

(Same  file  name  that  was  used  in  Line  130  in  file  TP2FIL.) 

Type  in 
RUN 

The  system  will  respond  with 
SNUMB  XXXXXT 


41.  After  the  finite  element  grid  data  have  been  prepared  according  to 
the  above  procedures,  the  grid  may  be  plotted  by  use  of  the  WES  Graphics  Com¬ 
patibility  System  (GCS)  (1979)  subroutines.  The  plot  file  can  be  directed 
either  to  a  Calcomp  plotter  or  to  a  Tektronix  4014  terminal  interactively. 

This  program  FNGRID  listing  is  presented  in  Figure  22. 

Steps  to  run  FNGRID 

42.  Program  FNGRID  may  be  run  by  using  the  following  procedures. 

A.  Use  the  XEDIT  mode  editor  to  change  the  total  number  of  nodal 
points  and  elements  of  the  grid. 

a.  Sign  on  CYBERNET  (Control  Data  Corporation)  COMPUTER  SYSTEM 

Telephone  No.  2047  for  300  Baud 
Telephone  No.  2030  for  1200  Baud 

b.  Call  out  the  procedure  file  PL0TF 
Type  in 

GET , PLOTF /UN=CER0M0 
SAVE,  PLOT F 

c.  Use  XEDIT  to  make  necessary  changes 
Type  in 

XEDIT 

d.  Type  in 
L/IXY/ 

The  system  will  respond  with 

Q , FNGRI D . C/ I XY/ 35 /* ; C/NEL/ 48/* 

e.  Change  the  total  number  of  nodes  (IXY)  and  elements  (NEL) 

to  the  desired  number 
Type  in 

C/35/TOTAL  NO.  OF  NODES/ 

C/ 48/TOTAL  NO.  OF  ELEMENTS/ 

f.  Option:  If  the  user  does  not  wish  to  change  data  file 

names  FINDAT  and  FINDAT1  to  other  names,  go  to 
section  g.  Otherwise, 

Type  in 

C/FINDAT/NEW  FILE  NAME/ 

C/FINDAT1 /NEW  FILE  NAME/ 

g.  Type  in 
E,  ,RL 

(This  command  will  terminate  the  editing  mode  and  the 
file  PLOTF  will  be  replaced  by  the  edited  version.) 
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OLD i FNGRID 
/LIST 

PROGRAM  GRIDf L< INPUT (OUTPUT. TAPES  -  INPUT. TAPF6. TAPES. 
4TAPE12.TAPE13.TAPE10.TAPE66.TAPE77.T APE9.TAPE99) 

C  ***  IXY  IS  THE  MAX  NUMBER  OF  NODES 
C  ***  MEL  IS  THE  MAX  HUMBER  OF  ELEMENTS 
C  *»*  N1-N4  IS  THE  NODE  NUMBER  OF  EACH  ELEMENT 
C  ***  DEP  IS  THE  DEPTH  AT  EACH  NODE 
C  *»*  10=0  PLOT  NODE  AMD  ELEMENT  NUMBERS 
C  ***  IC1  1  PLOT  ELEMENTS  NUMBERS. 

C  ***  101=2  NO  NUMBERS  PLOTTED 
C  ***  JC2-0  NO  DEPTH  DATA  OR'  PLOT 

C  ***  1C2=1  READ  AND  PLOT  DEPTH  DATA 

C  »■**  XTR.YTR-X  AND  V  COORDINATES  OF  THE  CFNTER  OF  THE  SF.HICIRCLE. 

C  *»»  ANGG-ANGLE  IN  RADIANS  OF  SEMICIRCLE  DIAMETER  FROM  GRID  X-AXIS. 

C  ***  RR  -RADIUS  OF  SEMICIRCLE  IN  FEET 

C  ***  NDEV.  0  PLOT  FILE  IS  DIRECTED  TO  THE  CALC0I1P  PLOTTER. 

C  ***  NDEV- 1  PLOT  FILE  IS  DIRECTED  TO  THE  TEKTRONIX  4014  TERMINAL. 

C  ***  NODR-NO •  OF  POINTS  OH  THE  CIRCUMFERENCE  OF  THE  SEMICIRCLE  (OR  CIRCLE). 

C  ***  NCL=  1  SEMICIRCLE  GRID  IS  USED. 

C  ***  NCL  0  FULL  CIRCLE  GRID  IS  USED. 

C  ***  FACT  IS  THE  FACTOR  USED  TO  INCREASE  OR  DECREASE  THE  ORIGIONAl.  PLOT  SIZE. 
C  ***  FACT- 1.0  IF  NDEV  1 

C  ***  UXMAX-  THE  R IGHT -MOST  WINDOW  BOUNDARY.  (SEE  GCS  MANUAL  ) 

C  ***  UXMIN-  THE  LEFT-MOST  WINDOW  BOUNDARY.  (SEE  GCS  MANUAL) 

C  ***  WYMAX-  THE  TOF-NOST  WINDOW  BOUNDARY.  (SEE  GCS  MANUAL) 

***  W Y M I N  =  THE  BOTTOM-MOST  WINDOW  BOUNDARY.  (SEE  GCS  MANUAL) 

***  XLO--  THE  X  COORDINATE  OF  THE  LOWER  LEFT  CONNER  OF  THE 
FIRST  CHARACTER  OF  THE  OUTPUT.  (SEE  DCS  MANUAL) 

***  YLO  THE  Y  COORDINATE  OF  THE  LOWER  LEFT  CONNER  OF  THE 
FIRST  CHARACTER  OF  THE  OUTPUT.  (SEE  GCS  MANUAL) 

***  RTITLE--  TITLE  OF  THE  PLOT 
***  NAME l-  IDENTIFIER  USED  FOR  THE  PROJECT 
***  USE  PROC.  FILE  *  P  L  0  T  F  *  TO  SET  DESIRED  DIMENSIONS. 


COMMON  LDUM ( 2000 ) 

DIMENSION  X ( IXY) iY(IXY) ,DEP(IXY) 

DIMENSION  K ( NEL ) . N 1 <  NEL ) i N2  <  NEL  >  >N3 ( NEL ) »N4(NEL) 
DIMENSION  XAVE ( MEL  > . VAUE(MEL) 

DIMENSION  N A M E 1 ( 10) 

DIMENSION  RT I TLE ( 8 ) 

R  E  A  D ( 12.10)  NDEV. NODR. NCL. FACT 
10  FORMAT (315. FS.  1  ) 

I R 1  -  •  R 1 ' 

ZVt  0 
IV2  -1 

CALL  SET JCI ( IR1 . I  VI  ) 

IF ( NDEV . EC . 1 )  CALL  SET JC I ( I R 1 . I V2 ) 

READ (12. 10)  NAME1 
15  FORMAT  < 1  OP 6 ) 

NSPACE  =N0DR-1 
READ< 12.20)  SCALE. IC1 . IC2 
20  FORMAT <  F 10 . 2 »  21 10  ) 

READ (12. AO)  XTR.YTR. ANGG.RR 
40  FORMAT ( 4F 10 . 2 ) 

READ (13.50)  (X( I ) . Y ( I  )  .  I  -1 .  IXY  ) 

50  FORMAT ( 3( 6X  >  2F3 . 3 ) ) 

FA  I  =3.141o 

IF (NCL. EC. 1)  DAO  PAI /NSPACE 
IF (NCL.EO.O)  D  A  2  =  2.«PAI/N0DR 
NNOD ' IXY 
N0DB2  NNOD-NODR 


Figure  22.  Program  FNGRID  for  plotting  the  finite  element  solution 
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U  a  i i £ i c i > ^ , 

UR  I T  E ( 6  > 1 05  >  NAME1 
URITEC8.96) 

95  FORMAT!//.  iOX*K«***»  4  *.M:***4:*******it* *■»•*•*****«#*****«****»**<*•  ) 
105  FORMAT ( / >  20X , 1 0 A6  >  7 ) 

96  FORMAT (10X»***«*<**«4 **#***«**«*»*»#*** *****«*#« ***F****i«*4U'»/) 

C  YY*  l\  (  I  )  USE  TEMPORARY 

DO  60  1  l.IXY 
60  K  < I  1  I 

F:  E  A  P  (  9  >  7  0  >  (PEP  Cl )  .  I  -l  .((NOD) 

70  FORMAT! 16F5. 2) 

URI TE ( 8 i 1 1 5 ) 

115  FORMAT (7X. • RODE  * , 1 IX, ' X.-COOP* . 11X. * Y-COOD' , 1 1 X . • DEPTH  * ) 

URITE(8,90>  ( X ( I > . X ! I ) >  Y ( I ) >  D  E  P ( 1 ) . I  l.IXY) 

90  F0RHAT(I10»2X,F15.2,2X,F15.2,2X,F15.2> 

NO  DP l -  IXY-NODR+l 

N  =  I  X  Y 
M  =  NEL 

DO  130  I  l.N 
X< I )  X ( I ) f SCALE 
Y ( I ) . Y( I >» SCALE 
130  CONTINUE 
C  *«i  PRINT  INTERMEDIATE  DATA 
1  AO  FORMAT ( IX. 12A6 ) 

URI TE(8r 150)N»M. YRr SCALE 

150  FORMAT ( IX. 12HH0  OF  POINTS. 15./. IX. 1  9HN0  OF  ELEMENTS . 15 ./. IX . 1 6HRAN 
*GE  OF  Y -COORD. F10.0./. 1X.22HSCALE  FACTOR  (IM/FEET ) . E 10 . 3 ) 

C  ***  CALCULATE  ELEMENT  CENTROIDS 
UR  I  TEC  8. 160) 

160  FORMAT (// . 7HELEMENT . 1 X .  -1HN ( 1 ) . 2X . 9HNC2) . 2X . 1HNC3) . 2X . 7HELEMENT , IX, 
*9HN( 1 ) ,2X, 9HNC2) .2X. 9HNC3) . 2X , 7HELEMENT , IX,  1HN( 1 ) , 2  X , 6  H  N ( 2 ) » 2X , 9HN 
*  (  3  )  ) 

READ C 9 , 170)  ( N1 ( I ) , N2 ( I ) , N3 ( I ) . I  - 1 , NEL  > 

170  F0RMATC3C5X.3I5)) 

ISO  FORMAT ( A 16, 2 X, 616, 2 X, 416) 

URITE (8, 130 ) < I »H1 < I ) .H2C I ) ,N3( I ) .1  1 , MEL ) 

URI TEC  10,250) CXC I) , I  -l.IXY) 

UR  I  TEC  10,250) CYC  I) , I  -l.IXY) 

UR  I  TEC  10, 170)  CK'l  ( I  >  ,N2C  I )  ,N3(  I )  .  I  =1 ,  NEL) 

UR  I  TEC  10, 250)  (DEPCI >  ,  I  -  1 , NMOD )  ' 

250  F0RMATC8F10.2) 

C  **»  READ  IN  DATA  FOR  X-RANGE  AND  Y-RANGE  *** 

READ  <  12,251  )  U  X  MAX  ,  U  XIII M  ,  UY  MAX  ,  UY  M  I N 

251  FORMAT ( 9F 1 0 . 1 ) 

C  ***  READ  IN  DATA  TO  LAPEL  THE  PLOT C  *** 

REAP (12, 61)  XLO.YLQ 

91  F0RMATC2F6.2) 

READ!  12,62)  RTITLE 
URITE(B,92>  RTITLE 

92  FORMAT  C8A10) 

DO  35  1=  1  ,M 

35  N9 (  I  )  0 

IF  C IC1 .CE. 2)  GO  TO  220 
DO  200  I  1,H 
NM1  NHI) 

NN2  M 2  C  I ) 

NN3’ N3< I ) 

NN9  F'  4  C  I  ) 

X  AVE  < I  >•  XCNN1 >  +  XCNN2)+XCNN3> 

YAVE  (  I  )  =  Y  (  HU  )  +Y  NN2  )  +  Y  (  NN3  > 

IF  <N6( I ) .GT .0)  GO  TO  210 
X  AVE  <  I  >  ):AVE  (  1  ,  /  3  . 

YAVE ( 1 >;  YAVE ( I ) / 3 . 

GO  TO  200 

210  X  A  V  E  <  1  )  =  <  XAVE  <  I  )  tUIIHI  >  ) /A  . 

Y  AVE  <  !  >  <  Y  AVE  (  1  )  •»  Y  (  HN  \  )  )  .  A  . 
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200  C  _n  T  1  riut. 

220  COLL  US  TORT 

CALL  UFSETCSPEED*.l20.  > 

CALL  U  E  ft  A  S  E 
CALL  UDELL 

:  CALL  UFSET ( 'TERMINATOR ‘ > ' i " ) 

IFUIDEV.FQ. 1 )  GO  TO  225 
CALL  USE  T  < *ALTE‘  ) 

CALL  UEKASE 
225  CONTINUE 

CALL  USE T  < •DEVICE*  ) 

CALL  FACTO ft (FACT) 

CALL  UDAft£A< 0 . 0. 14 . 0.0 .5 > 10 . ) 

CALL  USE T ( "VIRTUAL* ) 

CALL  UU I NDO ( WXMAX » UXMI N . MYMAX  >  UYMIN ) 
DO  25  1=1. h 
XI=  I 

C  »**  HOVE  FEN  TO  FIRST  ELEMENT  NODE 
Ml  =N1 < I > 

M2  -N2  (  I  ) 

M  3  -  M  3  (  I  ) 

M4MT4  (  I  ) 

CALL  USET ( 'NQAXES*  > 

CALL  USETCLIME*) 

CALL  UN0VE(X(H1  )  .  Y(M1 )  ) 

CALL  UFEN ( X ( M2 ) . Y  <  M2 ) ) 

CALL  UNO  VE < X  <  M2  >  >  Y  <  N2  > ) 

CALL  UF  EN ( X  <  M3 ) . Y ( M3 ) ) 

CALL  UHOVE ( X ( K3 ) »  Y ( M3) ) 

IFUH.EO.O)  GO  TO  26 
CALL  UFEN ( X(M4).Y(M4)) 

CALL  UMOVE ( X  <  N4 ) . Y ( M4 ) ) 

26  CALL  UPENC<<iU  )  >Y(M1  )  ) 

25  CONTINUE 

I F ( I C 1 . EG  .  2  )  GO  TO  27 
DO  275  I '  1  >  H 
ELE  I 

CALL  UHOVE ( XAVE ( I ) »  YAVE  < I ) ) 

CALL  USE  T ( *  SMALL  * ) 

CALL  USETt *  INTEGER*  ) 

CALL  UFftINT ( XAVE ( I ) , YAVEi I ) . ELE ) 

273  CONTINUE 

IF  <  IC1  .  GE  .  1  )  GO  TO  27 
DO  270  I  I  ,N 
EN=  I 

CALL  UHOVE ( X ( I ) i Y ( I ) > 

CALL  USE  T  <  ’SMALL* ) 

CALL  USE T  <  *  INTEGER’  ) 

CALL  UFftINT <X( I ) . Y ( I ) .EN> 

270  CONTINUE 

27  CONTINUE 

IF  < 102 , EC. 0)  GO  TO  116 
CALL  USET ( ‘HARDWARE’ > 

CALL  USET ( ’SMALL’ ) 

DO  117  I-l.NNOD 
CALL  USET ( ’REAL* ) 

CALL  Uf  R I  NT ( X  < I ) >  Y ( I ) >  DEF< I ) ) 

117  CONTINUE 
116  CONTINUE 
C  I*  *  L ASEL  FLO  T  *** 

CALL  USE  T  <  'VIRTUAL ’ ) 

CALL  USE  T  <  *  L  A  ft  G  f"  *  ) 

CALL  USET< ’TEXT* I 

CALL  UFftINT (XLO.YLO.RTITLE) 

CALL  VEND 
UF.  1  T  F  (  A  »  l  0  00  ) 

1000  FORMAT  (  IX.  *F  LOTTING  LOl.FLEIt*) 

STOP 

END 
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Use  procedure  file  to  run  program  FNGRID. 

a.  Type  in 
OLD , PLOTF 
BEGIN, RUNPL,PLOTF 

The  system  will  respond  by  printing  out  the  statements 
which  are  edited  by  "Q"  statement  in  file  PLOTF 

b.  Answer  the  first  question  of  PROGRAM  FILE  NAME 
Type  in 

FNGRID1 

c.  Answer  the  question  of  TYPES  (X,N,B)  OF  FILE 
Type  in 

S 

d.  Answer  the  question  of  DEVICE 
Type  in 

DR4  (if  Calcomp  plotter  is  chosen) 

TK4  (if  Tektronix  4014  terminal  is  chosen) 

(Note:  Computer  will  begin  to  plot  on  the  screen  if 

TK4  is  chosen;  computer  will  respond  with  more  ques¬ 
tions  if  DR4  is  chosen.) 

e.  Answer  the  question  of  USER'S  ID 
Type  in 

CEROXX  (the  user's  ID) 

The  system  will  respond  with 
TAPE  99  ROUTED  THRU  AJZZ123 

(where  AJZZ123  is  the  job  name  of  the  plot  file) 

f.  Call  (Ext  3442)  the  operator  of  the  COPE  terminal  at  the 
WES  Automatic  Data  Processing  (ADP)  Center  to  retrieve  the 
plot  file.  (Note:  TAPE10  is  saved  and  replaced  by  file 
name  FINDAT  after  the  run,  and  will  be  used  as  data  file  to 
run  program  FINITE.) 


Data  preparation  information 


follows . 


The  variable  names  and  format  allocations  for  program  FNGRID  are  as 


DATA  SET  NO.  1 


CARD  NO.  1 


Format  Variable  Name 


1  -  5 


6-10 
11  -  15 


_ Description _ _ 

NDEV=0 ,  plot  file  is  directed  to  the  Calcomp 
plotter 

NDEV=1,  plot  file  is  directed  to  the  Tek¬ 
tronix  4014  terminal 

Number  of  points  on  circumference 

NCL=0,  full  circle  grid  is  used 
NCL=1,  semicircle  grid  is  used 


NCL 


Col. 


Format  Variable  Name 


Description 


16  -  20 

F5  .  1 

FACT 

Factor  used  to  change  original  plot  size. 

(Note:  If  NDEV-1,  then  FACT=1.0.) 

1  -  60 

10A6 

NAME1 

CARD  NO.  2 

Identifier  used  for  the  project 

~  • 

■ 

CARD  NO.  3 

.  .  «  ,  ■ 

1 

1  -  10 

F10.2 

SCALE 

SCALE=1  and  2,  x  and  y  coordinates  of 
the  nodes  are  in  feet  and  inches, 
respectively 

• 

11-20 

no 

IC1 

IC1=0,  plot  node  numbers  and  element  numbers 
IC1=1,  plot  element  numbers  only 

IC1=2,  no  number  plotted 

t 

21  -  30 

no 

IC2 

IC2=0,  no  depth  data  on  plot 

IC2=1 ,  read  and  plot  depth  data 

• 

CARD  NO.  4 

1  -  10 

F10.2 

XTR 

X  -  coordinate  of  the  center  of  semicircle 

!» 

11-20 

F10.2 

YTR 

Y  -  coordinate  of  the  center  of  semicircle 

.  • 

21-30 

F10.2 

ANGG 

Angle  of  the  semicircle  (radians) 

31  -  40 

F10.2 

RR 

Radius  of  the  semicircle  (feet) 

CARD  NO.  5 


1  -  10 

F10 . 1 

WXMAX 

11-20 

F10 . 1 

WXMIN 

21  -  30 

F10 . 1 

WYMAX 

31  -  40 

F10 . 1 

WYMIN 

1  -  6 

F6.2 

XLO 

7-12 

F6.2 

YLO 

The  rightmost  window  boundary  (See  GCS  1979 
manual) 

The  leftmost  window  boundary  (See  GCS  1979 
manual) 

The  topmost  window  boundary  (See  GCS  1979 
manual ) 

The  bottommost  window  boundary  (See  GCS  1979 
manual ) 

CARD  NO.  6 

x-coordinate  of  the  lower  left  corner  of  the 
first  character  of  the  output  (See  GCS 
1979  manual) 

y-coordinate  of  the  lower  left  corner  of  the 
first  character  of  the  output  (See  GCS 
1979  manual) 


Col. 


Format  Variable  Name 


Description 


CARD  NO.  7 


1  -  80 

8A10 

RT1TLE 

Title  of  plot 

DATA  SET  NO.  2 

CARD  NO.  1 

(Note : 

This  data 

set  is  obtained  after  the  program  CONVER  is  run.) 

1  -  6 

Blanks 

7-14 

F8.3 

X(l) 

x-coordinate  of  node  1 

15  -  22 

F8.3 

Y(l) 

y-coordinate  of  node  1 

23  -  28 

Blanks 

29  -  36 

F8.3 

X(2) 

x-coordinate  of  node  2 

37  -  44 

F8.3 

Y(2) 

y-coordinate  of  node  2 

• 

45  -  50 

Blanks 

51  -  58 

F8.3 

X(3) 

x-coordinate  of  node  3 

59  -  66 

F8.3 

Y(3) 

y-coordinate  of  node  3 

(Note : 

Six  values  on  each  card 

Total  number  of  x  and  y  values  should  equal 

m 

the  tota 

1  number 

of  nodes.) 

DATA  SET  NO.  3 

CARD  NO.  1 

1  -  5 

F5.1 

DEP(l) 

Depth  in  feet  at  each  nodal  point 

v~ 

6-10 

F5 . 1 

DEP(2) 

(Note:  Sixteen  values  on  each  card,  and  use 

as  many  cards  as  necessary.) 

76 

-  80 

F5 . 1 

DEP( 16) 

DATA  SET  NO.  4 

• 

CARD  NO.  1 

1 

-  5 

15 

K(  1 ) 

Dummy  variable  used  as  an  index  number  for 

the  element 

6 

-  10 

15 

Nl(l) 

First  nodal  point  of  element  1 

11 

-  15 

15 

N2  ( 1 ) 

Second  nodal  point  of  element  1 

16 

-  20 

15 

N3(l) 

Third  nodal  point  of  element  1 

-“.V/ 

39 


Col. 

Format 

Variable  Name 

Description 

21  -  25 

15 

K(2) 

Dummy  variable 

26  -  30 

15 

Nl  (2) 

First  nodal  point  of  element  2 

31  -  35 

15 

N2(2) 

Second  nodal  point  of  element  2 

36  -  40 

15 

N3(2) 

Third  nodal  point  of  element  3 

41  -  45 

15 

E(3) 

Dummy  variable 

46  -  50 

15 

Nl  (3) 

First  nodal  point  of  element  3 

51  -  55 

15 

N2(3) 

Second  nodal  point  of  element  3 

56  -  40 

15 

N3(3) 

Third  nodal  point  of  element  3 

(Note : 
Nl,  N2 , 

Twelve 
and  N3 

values  on  each  card  for  as  many  cards  as  needed.  The  order  of 
of  each  element  is  chosen  by  following  a  counterclockwise 

direction. ) 

Example  run 

44.  The  finite  element  grid  of  Figure  17  will  be  processed  as  a  typical 

representative  example  of  the  capability  of  the  programs. 

A.  Use  XEDIT  mode  to  set  the  total  number  of  nodes  to  35  and  total 
number  of  elements  to  48  (Figure  23). 

GET » PL0TF/UN=CER0M0 
/SAVE » PLOTF 
/OLD » PLOTF 
/XEDIT 

XEDIT  3.1.00 
??  L/Q » FNGRID/ 

G» FNGRID. C/I XY/30/* 5 C/NEL/40/* 

??  C/30/35/ 

Q»FNGRIH.C/IXY/35/*? C/NEL/40/* 

??  C/40/48/ 

0»FNGRID.C/IXY/35/*»C/NEL/4R/* 

??  E»»RL 

PLOTF  REPLACED 
PLOTF  IS  A  LOCAL  FILE 
AEB  »  0.234UN7S. 

/ 

Figure  23.  Using  XEDIT  mode  to  set 
total  number  of  nodes  and  elements 


B.  Run  program  FNGRID 

a.  Plot  file  is  directed  to  the  Tektronix  4014  terminal 
(Figure  24) . 


40 


OLD.PLOTF 

/ICCtN.RUNPL.PLOTF.FDl  .FD2.FD3. 

C  III  35  IS  THE  MAX  NUM1ER  OF  NODES 
DIMENSION  XOS  >,  V(  3S  J,D£P(35 ) 

READ! 13.58 >  (X<l).Y<I).I*1.3S> 

NHOD-3S 
DO  6*  1*1,35 

WRITE!  8.90)  (KII).XII),V!I ),  DEPCI  ),I*1,3S> 

N0DII-3S-N0DR* 1 

H*3S 

WRITE!  18,254  >(X(I), 1*1,35) 

WRITE!  10.258  XV!  I),  1*1, 3S> 

END  OF  FILE 

C  111  48  IS  THE  RAX  NUN8ER  OF  ELEMENTS 

DIMENSION  K ( 48  ) , N1 1 48  > ,N2( 48  >,N3(4B  >,N4<48 > 
DIMENSION  XMC(48).VAUC!48) 

M*48 

READ! 9, 178 >  INI !I ),N2! I >,N3!I ), 1*1,48  ) 

URITE(8.18S)II,N1II),H2II).N3II),I*1,48) 

URITEI18,178)!N1CI),N2II>,N3II),I*1,48) 

END  Of  FILE 

ENTER  PROGRAM  FILENAME 
7  FNCRID1 

ENTER  S  FOR  SOURCE  FILE  WITH  NO  LINE  NUMBERS 
ENTER  N  FOR  SOURCE  FILE  WITH  LHC  NUMIERS 
ENTER  I  FOR  IINARV  FILE 
?  S 

DEUICE- 

?  TK4 


Figure  24.  Directing  plot  file  to 
Tektronix  4014  terminal 


b.  Plot  file  is  directed  to  the  Calcomp  plotter  (Figure  25), 

OLD.PLOTF 

/BEGIN, RUNPL.PL0TF.FDl.FDa.FD3. 

C  **t  35  IS  THE  flAX  NUMBER  OF  NODES 
DIMENSION  X(35  ),V!35  >,DEP(35  ) 

READ! 13, S9 )  (XI I ), V! I ), I • 1 , 35 ) 

NN0D-3S 
DO  60  1*1.35 

URITEI8.90)  (K!I),X(I>,V(l),DEPm.I*1.35I 

NODB1-3S-NODR41 

N*35 

WRITE! 10,259  )(X(I),I*1,35) 

WRITE! 10,250  )!Y!I),I*1,35> 

END  OF  FILE 

C  8*1  48  IS  THE  MAX  NUMBER  OF  ELEMENTS 

DIMENSION  K(48),N1!48)>N2(48),N3(48),N4(48) 

DIMENSION  XAUE148  ), VAUEI48 ) 

M-48 

READ(9. 170 )  <N1(I),N3(I),N3(I),I*1,48) 

WRITE(8,180)(I,NHI  ),N2!  I  >,N3!  I ),  1*1,48) 
UR1TE!10,170)!N1!I),N2!I),N3!I),I*1.48) 

END  OF  FILE 

ENTER  PROGRAM  FILENAME 
7  FNCRID1 

ENTER  S  FOR  SOURCE  FILE  WITH  NO  LINE  NUMBERS 
ENTER  N  FOR  SOURCE  FILE  WITH  LINE  NUMBERS 
ENTER  B  FOR  BINARY  FILE 
7  5 

DEUICE- 
7  DR4 

ENTER  YOUR  USER  NUMBER 
7  CER0N8 

TAPE99  ROUTED  THRU  JOB  AJZZ324 
CALL  DOT  GRIFFIN  < TEL. 3442)  k  SPECIFY  3-BIT  SOFTWARE 
REUERT. 

/ 

Figure  25.  Directing  plot  file  to  Calcomp 
plotter  terminal 


C.  Listing  of  data  files  (FD1,  FD2 ,  FD3)  (Figure  26). 


OLD » FD1 
/LIST 

0  9 

1  1.0 

SAMPLE  PROBLEM 

1.00 

0  0 

0.00 

0.00  0.00 

-35.0 

35.0  -40.0 

. -8.00-37.00 

FINITE  ELEMENT 

GRIDS 

Figure  26a. 

Listing  of  data  file  FD1 

0LD»FD2 

/LIST 

-8.000  -32.000 

0.000  -32,000 

-4.000  -28.000 

4.000  -28.000 

0.000  -24.000 

8.000  -24.000 

4.000  -20.000 

-8.000  -16.000 

8.000  -16.000 

-4.000  -12.000 

-8.000  -8.000 

0.000  -8.000 

-4.000  -4.000 

4.000  -4.000 

0.000  0.000 

8.000  0.000 

4.000  4.000 

0.000  8.000 

11.000  4.600 

8.500  8.500 

0.000  12.000 

-4.600  11.100 

-11.100  4.600 

-12.000  0.000 

9.00 

15.0 


8.000  -32.000 
-8,000  -24 . 000 
-4.000  -20.000 


0.000  -16.000 
4.000  -12.000 
8.000  -8.000 
-8.000  0.000 
-4.000  4.000 

12.000  0.000 

4.600  11.100 

-8.500  8.500 


Figure  26b.  Listing  of  data  file  FD2 


0LD.FD3 


/LIST 

0.20 

0.20 

0.20 

0.20  0.20  0. 

20 

0.20 

0.20  0.20  0.20 

0.20 

0. 

2'.j 

0.20 

0.20 

0.20 

0.20  0.20  0. 

20 

0.20 

0.20  0.20  0.20 

0.20 

0. 

20 

0.23 

0.23 

0.25 

I 

2 

4 

2 

3 

3  1 

4 

6 

4 

7 

6 

n 

7 

4  2 

5 

7 

5 

8 

7 

3 

8 

3  6 

7 

9 

7 

8 

10 

6 

9 

11  9 

12 

11 

7 

12 

9 

7 

10 

12  10 

13 

12 

8 

13 

10 

11 

12 

14  12 

13 

10 

11 

14 

16 

14 

17 

16  12 

1  7 

14 

12 

IS 

17 

15 

18 

17  13 

IB 

15 

16 

17 

19 

17 

18 

20  16 

19 

21 

19 

21 

17 

22 

19  17 

20 

*1  *1 

20 

23 

22 

IB 

23 

20  35 

21 

34 

21 

22 

24 

*1  n 

23 

20  23 

27 

28 

21 

24 

34 

22 

23 

24  23 

28 

25 

34 

24 

33 

24 

23 

26  25 

28 

29 

24 

32 

33 

24 

26 

32  23 

30 

26 

25 

29 

30 

26 

31 

32  26 

30 

31 

Figure  26c. 

Listing  of  data  file  FD3 

0.2: 


42 


The  computer  plot  of  the  finite  element  grid  originally  conceived  in  Figure 
is  presented  in  Figure  27  after  processing  by  the  programs  DIG],  TP2FIL, 
CONVER,  and  FNGRID . 


Figure  27.  Computer  plot  of 
finite  element  grid 


PART  IV:  PROGRAM  FINITE 


45.  The  propagation  of  periodic,  small  amplitude  surface  gravity  waves 
over  a  variable  depth  seabed  of  mild  slope  is  governed  by 


A  "  (cc  Ad)  +  l2d  =  0 

g  c 


(1  bis) 


Equation  1  is  solved  numerically  by  program  FINITE  using  a  hybrid  finite  ele¬ 
ment  method  originally  developed  by  Chen  and  Mei  (1974)  to  solve  the  diffrac¬ 
tion  Helmholtz  equation  in  a  constant  depth  region. 


Data  Preparation  for  Program  FINITE 

46.  Two  types  of  data  are  necessary  to  run  program  FINITE:  (a)  data 
values  defined  by  parameter  statements  in  the  program  FINITE  and  (b)  data  sub¬ 
mitted  to  the  program  FINITE  by  data  files.  Program  FINITE  is  presented  in 
Figure  28. 

Data  values  defined 
by  parameter  statements 

47.  The  following  data  values  are  defined  by  parameter  statements  in 


the  program  FINITE. 

Variable  Name 

Format 

Description 

NNOD 

Integer 

Total  number  of  nodes 

NELE 

Integer 

Total  number  of  elements 

NODE 

Integer 

Total  number  of  radiation  boundary  nodal  points 

NCS- 

Integer 

Total  number  of  coefficients  in  the  expansion  of 
radiation  domain 

ISPT 

Integer 

Number  of  selected  elevation  station  (nodes) 

IS  PTE 

Integer 

Number  of  selected  elevation  station  (elements) 

NSYSK 

Integer 

NLGL  x  NBAND  x  2  +  50 

(see  comment  statements  in  program  FINITE) 

NSYTP 

Integer 

NLGL  x  NBAND  +  100 

(see  comment  statements  in  program  FINITE) 

('•'■Note:  The  value  of  NCS  is  not  known  a  priori.  Enough  terms  must  be  included 
such  that  additional  terms  have  a  negligible  effect  on  the  solution.  NCS  can 
be  estimated  by  determining  the  argument  (2^/0  of  the  J  Bessel  Functions 


(Committee  for  the  Calculation  of  Mathematical  Tables  1958).  The  values  of 
different  order  J  Bessel  Functions  with  this  argument  can  be  determined  using 
Bessel  Function  tables.  When  the  value  of  a  Bessel  Function  of  higher  order 
becomes  much  less  than  the  lower  order  Bessel  Function  of  the  same  argument, 
NCS  is  set  equal  to  the  order  of  the  higher  order  Bessel  Function.  For 
h  =  1  ,  k  =  12  ,  T  =  2  sec;  1  =  (Vgh)T  =  11.4  ft,  (2rtR/L)  =  6.64  .  Since 
Ji3^6'64)  is  much  less  than  Jq(6.64)  ,  NCS  is  set  equal  to  13.  The  model 
should  then  be  run  with  NCS  equal  to  both  13  and  14.  If  the  results  are  the 
same  within  desired  tolerances,  the  value  of  NCS  =  13  is  confirmed.) 

Data  values  submitted  by  data  files 

48.  Data  file  FINDAT1  contains  Data  Set  No.  1  through  Data  Set  No.  3. 
Data  file  FINDAT  contains  Data  Set  No.  4  through  Data  Set  No.  7.  (FINDAT  can 
be  obtained  after  running  the  plotting  program  FNGRID.) 


OLD* FINITE 
/LIST 

*DECK  MAIN 

P  ROGRAM  L ALB  <  I NPUT  .  OUTP  UT ,  T APE5  =  I NFUT  ,  TAPE6  =OUTPUT  r  T APE4  , 
1TAFE9.TAFE11  ,  T  APE  1 2 .  T  APE  1 3 .  T  APE  1 4  ,  T  APE1S .  T APE  1 6 . T APE  17 . T AFE 1 8 , 
1TAPE19.TAPE20.TAFE21.  T AF  E22 ,  T APE23 ,  TAPE2  4  »  TAPF.23  > 

C 

CCC  BARBERS  POINT  HARBOR  OSCILLATION  STUDY 
C 

C  ******  »:*»***  **  **************************  ****  *********************** 

C  MELE  '  TOTAL  NO.  OF  ELENEMTS 

C  NNOD*  TOTAL  NO.  OF  NODES 

C  NODR  TOTAL  MO. OF  R AD  I  AT  I  ON-BOUNDARY  NODAL  POINTS 

C  NCS  -TOTAL  NO. OF  COEFF.  IN  THE  EXPANSION  OF  RADIATION  DOMAIN 

C  HEQT-  total  MO.  OF  SIM.  LIN.  EOS.  = NHODFNODR 

C************»******************************************************** 
COMPLEX  SYSK.SYSO.XH.DH.SM 
*»  SYTP . STKP » FTEHP 

COMMON/DB/  UK. ALPHA, D1 .DA2.DA4.RR, UR. XLiUl.  , RKHP.HKL 
DIMENSION  X ( 1277 ), Y ( 1277 ), NOD ( 1277 ), HCOM ( 2334 » 3  >«  XJ ( 25 ) » 

1  XY(25)  .XH(25>  ,DH<25>  r  ANGL(13)  ,SM(13>  » HD (2334)  , HSP( 100 > , IPER ( 5 > 
DIMENSION  AF<  12?/)  .AT<  1297)  >  ID < 20  ) , NSF'E  (  100) 

DIMENSION  U<  2334 ) , V ( 2334 >, YC ( 2334 > ,Z<3> 

COHMON/Ll/SYSi<(  19850) 

C0MM0N/L2/SYS0< 1300) 

COKM0N/L3/STMP(130O> 

C0MM0M/L1 /SYTP (9950) 

DIMENSION  FTEMP ( 9950  ) 

DIMENSION  V J ( 105 ) , VY ( 105 ) 

DIMENSION  NR ( 3 ) 

DIMENSION  XK  < 25 ) 

COMMON  /NEU/  NUMBLK.NLG.NLGL.L11.LS 
1  .L1.L3 

LEVEL  2* SYSN.SYSOrSTHP. SYTP .FTEHP 
EQUIVALENCE  (FTEHP ( 1 ) . SYTP ( 1 ) ) 

DATA  ( I  PER ( J ) . J  - 1 . 5 ) /4HWAVE . 4H  PER . 4HI0D  .4HIH  S.4HECS  =  / 

I  PC  ••  0 
IF  MAX  =•  0 
R  E  A  D ( 5 . 2 )  ID 

2  FORMAT (20 A 4) 

WRITE (6. 3)  ID 

3  FOR M AT (//.35X. 20 A 4 >////) 

RE AD ( 5 » 5430 )NBAND»  NNOD . NELE . NODR. NCS. I  PR  » RR . ALPHA, DM . IBLKLEN 

5430  FORMAT ( 615 . F10. 2, F 10. 5. F 10. 2. 116) 

NELE3- 39NELE 

NEOT  -MNOD  +  NCS 
H0DR1  NODR-1 

NUMBLK  =  NEOT <  NBAND/ 1 BLKLEN  +  1 

PAI  -  3.141592654 

D 1  '  1.0 

DA2  -  FA  I /NODR  1 

DA4  =  DA2/2. 

5431  NLG-NEOT/NUMBLK 

NLGl.  NEOT- (NUMBLK-1 )*NLG 
nlgg  2*NLGL 

IF  (NLGGlNBAiKD.LT  .  198501G0  TO  5432 
NUMBLK  -  NUMBLK  t  1 
GO  TO  5431 

5432  CONTINUE 

URITE(6.  12)NELE.  I’NOD , HODR , NCS . NEOT , ALPHA , DA2 , DA4 , RR , D1 , DM 
12  FORMAT ( /IX . 23HT0TAL  NO.  OF  ELEMENTS  • » 15 « 3X , 20HT0TAL  NO.  OF  NODES 
1*.I5,3X»22HT0TAL  HO.  OF  SOURCES  =  , 1 5 , 3X , 30HF I RST  ELEMENT  W.R.T.  SO 
2URCES  I5./1X.30HT0TAL  NO.  OF  SIM.  LINEAR  EOS  * . 15. 4X . 25HUAVE  ANG 
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3LE  OF  INCIDENCE  -,F10.5,AX,29HINCHEMtNl  IS  AN..  ...KIUION  -  2K10. 
4,  /  iy.t  32rtOUTEST  EXTEND  Of  FINITE  DOMAIN  f  1  0  .  2  1  ;»hnOKM,  DFF  TM  0 
3F  OUTER  REGION  ■  ,  F3 . 2 , 5X . 3 1 HAVER .  DEFTH!FT>  OP  SmlCIRCLE  ‘.F10.2, 

6//  ) 

c<*»mimuiituim«i  HtittiiitmtmtmutTttmttitttmmmt 
UK  WAVE  HO. 

ALPHA  WAVE  ANGLE  Of  INCIDENCE 
01  WATER  PERTH  Of  OUTER  REGION 
DAO  INCREMENT  IN  ANGULAR  DIRECTION 
tltltllMUUtMMitllltlimHimiillMIIIIIIIIMIIIUmMMU 

IF  (IFF:  .ME.  1>  GO  TO  2013 
PRINT  2010 

2010  FORMAT k30X , 28Hli,F UT  DATA  CHEO-  EChO  M- I ■  ' 

FRUIT  201A,  NUMBLK  >MLGi  riLGL  •  I  ti.  RLEm  .  IP  R  .  nR.<n( 

201  -1  FORMAT!//. IX. 616) 

2013  CONTINUE 

REAp(5,300>TMAX, THIN. TDElT . IRR 
IF< IFF .EQ.OIGO  TO  3020 
PRINT  300.  TmA>:  .TMIN  .  TDELT  .  IKK 
300  FORMAT (3F1 0.2. 13) 

3020  XK(i)  =  TMIN 

DO  320  J  -  2. IRK 
J1  -  J  -  1 

XN ( J )  -  XK(Jl)  +  TDELT 

320  CONTINUE 

KEAD(3,305> IFC.IPMAX. ISF  T , IDPC. ISF  T£ , It  CH. 1BUGA. IBUGG. 1BUGL 
303  FORMAT (913) 

I  f  (  IF'R  .  EQ  .  0  )  GO  TO  3010 

PRINT  305.  IPC. IF' MAX. ISPT . ID PC . ISF TE. I  DC  Hi  1  PUG A . IBUGG. IDUGL 
3010  IF<IFC. EQ.OIGO  TO  3000 

RE AD <3. 306) <NSP<  J) .J'l.ISFT) 

RE  AD  (3. 306)  (NSPE<  J)  .  J-^l » ISPTE) 

306  FORMAT ( 1615) 

IF(IFR.EQ.O) GO  TO  3000 
URITE(6,306)  (NSF'(J)  .J-  l.ISPT) 

WRITE (6. 306) <NSPE< J) . 1 . 1SFTE) 

3000  IKT  •  IKK 

DO  99  IK  *  l.IKK 

WK  <2.»PAI)/<XK<IK)*SQRT<DH*32.2)) 

UR  I TE ! 6 «  73 )  :<K(IK),UK 

73  FORMAT (1H1 » // » 23X » 24HUAVE  PERIOD  IN  SECONDS  . FlO . 2 . 24Y. . 9HWAVE  NO 
1* » EI2 . 5  > // ) 

IF  <  IK  .GT.  1)  GO  TO  98 
CALL  I NPUT2  <  X.Y.NOD.NCON 

1  , N ELE.NNOD.NO DR, NCS.NEQT.N6AND . NCS2. DD . IDF'C.DM.IDCH) 

CALL  BAUD  (  .'ICON 

1 .NELE.NNOD. NODR, NCS.NECT.NB AND. NCS2) 

98  CONTINUE 

CALL  ASEMBK ( X , Y , NCOH , SYSK . SYSQ , S YTP , NLGL 
1,NELE.NN0D,N0DR,NCS,NE0T,NBAND.NCS2,DD,  NELE3 

I  .NLG.IDFC.IBUGA) 

IF(IBUGA.EO.O)GO  TO  3030 
WRITE <6, 330) 

350  FORMAT (//.1X.13HELK  ASSEMBLED) 

3030  CONTINUE 

CALL  GSPK ( SYTP » X J i XY  » XH .  PH » ANGL 
1  , NELE.NNOD.NODR.K'CS , NEOT . NB AND , NCS2 . NODR  1 , VJ, VY.NLCL, IBUGG) 

IF! IBUGG. EQ.0)G0  TO  3060 
WRI T  E ( 6 , 355 ) 

333  FORMAT.//, IX. 13HGSPK  COMPLETE) 

3060  CONTINUE 

CALL  LOAD ( SM , SYSQ , X J » DH • ANGL 
1, NELE.NNOD, NODR, NCS* NEOT. NBAND , HCS2. N0DR1 , IBUGL) 

IF< IBUGL. EQ.OIGO  TO  3070 
UR  I TE ( 6  »  373 ) 

373  FORMAT!//. IX. 13HL0AD  COMPLETE) 
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3070  CONTINUE 

L  ALL  EAHSOL  (  SY  $R  .  SVSC  .  STHF'i  NLG  >  >  K'UMBl.K  »  KEQT  » 

1  NLGG.SYTF  .  NLGL  .  F  TEMP  ) 

WF.  ITE  <  A.  20) 

20  F  ORMAT (////.ASX, 40<lH0)/»4OXf 

*40H  THE  SOLUTION  OF  THE  SYSTEM  /.ASX.AOCIHO)/// 

1  / 1 X  i  7HELEHEHT  >  OX  i  1  3H  AliF’L  I F  I C  AT  I  ON  ,  1 9X  .  8HUE  LOC I T  Y  1 1  4X  i  7HELEHENT  .  5X . 
213HAI1FLIFICATI0H,  1  9X  ,  8HUEL0C  I  T  Y  .  /  AX  ,  2H0R  .  1 1 X  .  6HF  ACT  OR  .  1 &)<  i 
3 1 OHCOMF  ONER'TS .  SX  » 9HH  AGHI TUOE  >  SX  i  2H0R .  1 1 X  i  6HF  ACTOR .  1 6X  , 

4  1  OHCOnF  QNE NT  S  >  8X  .  9HM AGN I  TUOE  >  /  2X »  AHNDE'E  .  AX  .  4 HM AG .  ,  BX  i  SHFHASE  »  9X » 

01  HU.  1 1  X  ,  1HV » 18X »  AHNOt'E  .  4X .  AHliftG .  ,  8X , SHF HASE , 9X , 1HU . 1 IX , 1 HV , /// ) 

Nl  -  HNOD+1 
00  22  1  1 .NEOT 
AR:  REAL  <  SYSO  <  I  >  ) 

A I  A I  MAG ( SY  SC (  I  )  ) 

AT  (  I  )  AT  AM2 ( A I , AR ) 

AF  ( I ) •  SGRT< AR*»2  +  AI**2) 

AF  <  I  >  AF  (  I  )  /  2  . 

22  CONTINUE 
C- 32 . 2 

00  A  0 1 4  L  l.NELE 
00  4  0  08  J-  1 . 3 
nF.  <  J  )  -  MOON  (  L  i  J  ) 

4008  CONTINUE 
1 1 -MR ( 1 ) 

I  2 -HR (2) 

1  3  -  N  R  (  3  ) 

XI  =  :<  (  I  1  ) 

X  2  -  X ( 12) 

X  3  - X ( 13) 

Y  1  •  Y  (  I  1  > 

Y  2  -  Y  l  1 2  > 

Y  3  -  Y  (  I  3  > 

H  -  Y3-Y2 
82-  Y1-Y3 

13  Y2-Y1 

Cl-  X  2 -  X  3 
02 *  X3-XI 
C3-  X 1  -  X  2 

AREA  0.G*<tl»C2-82*Cl) 

F  AC r  2 . YAREA/G 
L Uh  -  L 

LUM1  3TLUM-2 

LUH2  :  3  Y  L  U  M  ~  1 

L  U M 3  3YLUM 

PU1L1  1  -0 1  /F  AC 

0  U 1 L  2  --82/FAC 

0U1L3  1  -03/FAC 

0U2L1  * 'Cl/FAC 

0U2L2  --C2/FAC 

0U2L3  * -C3/FAC 

I F  (  AREA  . GT .  0.  )  GO  TO  4009 

URITE'.  4*4010)  AREA  .  LUM 

4010  FORMAT  (  IX.  21H0EKUG  ELI1K  100.  ARE  A  ‘ .  E  1 2 . 0 . 6H  AT.IS.UH-TH  ELEFIEN 

IT  ) 

STOP’ 

4009  COM  T I NUE 

SYTF(LUr,l>“  <0.0. 1 .0)*0V2L1 
S  Y  T  F  <  L  U  M 1 )  :  SYTFCLUMl)  +  DU1L1 
S  Y  T  F*  <  L  U  li  2  )  -'  <0.0.1  . O)»0V2L2 

SYTF<LUIi2)-  SYTF'LUK.2)  +  0  V  1 L  2 
SYTF  <LUit3)-  <  0 .0. 1 .0  )  Y0V2L3 
f,YTF(LUM3>-  SYTF<LUM3)  +  0V1L3 
4014  CONTINUE 
KU  0 

SIGMA-  2  .  *  F'AI/XM  IK  ) 
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DO  250  J=1 ,NELE 

HO  260  L=l,3 

N=NCON< J,L) 

60  2  (  L  >  •  AF  C  N) F2 . 

UC J>=0.0 
VC J>=0.0 
DO  270  L  =  1  >  3 
Kll-Kll+l 

DV1  ft  E  A  L  (  SYTFCKU)  ) 

DV2-  A I  MAG  (  SVTft(Ml)  ) 

U(J>  =  U(  J)+DV1»Z(L> 

70  V(J>-  VC J>+DV2»ZCL> 

UC  J)--UC  J)/Stl,MA 
VC  J)  -VC J)/S1GMA 

250  YCC J)=SGftT(U< JMUCJJ+VC J)*V< J) ) 

IFCIF'C.EC.OIGO  TO  502 
*WF\ITE<6,30) 

30  FORMAT  </>36X»4lHEL£VATI0N  AND  VELOCITY  FOR  SELECTED  ELEVATION  STAT 
1 1  OHS ( MODES ) .///) 

DO  500  1=  1  ,  ISf'T.2 
II  T  +  l 
J  =-  NSftCI) 

I F  <  1 1  .  GT  .  ISF'T )  GO  TO  501 
J1  -  NSFCJI) 

WRITE  (6.25)  J.AFC J) ,ATC J) ,UC J) »VC J) » YCC J) , J1 , AFC Jl> , ATC Jl) ,UC Jl>  * 
1VC  JD.YCC.il) 

GO  TO  500 

501  WRITE  <6. 24)  J  .  tiF  (  J  )  r  A  T  C  J  )  .  U  (  J  >  »  V  (  J  )  r  VC  (  J  ) 

500  CONTINUE 

502  NSK2  =  2YHSKIR 
IFCIREDRT.EQ.l >URI TE < 6 , 505 ) NSKIft 

505  FORMAT C /// , 39X  >  32HELEVATI0N  AND  VELOCITY  AT  EVERY  .I3.19H-TH  NODE 
10ft  ELEMENT,///) 

DO  400  I • 1 .NN0D.NSK2 
II  -*  I+NSKIft 
IFCII.GT.NN0D)G0  TO  410 

WRITE  C 4,25)  I  ,  AF ( I ) , AT ( I ) , U ( I ) > V < I ) , YC ( I ) , 1 1 , AF C 1 1 ) , AT C 1 1 > , U C 1 1 ) , 
1 V  <  1 1  >  ,  Y  C  <  1 1  ) 

GO  TO  400 

410  WRITE  <6.24)  I . AF < I ) , AT < I ) , U C I ) • V ( I ) , YC ( I > 

400  CONTINUE 

NNOI'l  NNOD+1 
UR  I TE  <  6 , 23 ) 

DO  414  I - NN0D1 . NEQT , 2 
IIs  I  +  1 

I F ( 1 1 , DT . NEQT  >  GO  TO  415 

UNITE  <6,25)  I  .AFC  I ) , ATC I ) ,U< I ) , VC  I ) , YCC  I  )  r II , AFC  II ) , ATC II ) ,UC II > , 
1VC II ) .YCCII > 

GO  TO  414 

4  15  WRITE  <6,24)  I . AFC  I ) , ATC I )  ,U< I ) »V< I ) , YCC I  ) 

414  CONTINUE 

UR  I TE  <  6 , 26  >  NSK I  ft 

26  F0RMATC////.43X.20H VELOCITIES  FOR  EVERY, I3.22H-TH  REMAINING  ELEMEN 

IT , ///  ) 

NECT1  'NEQT  + 1 
N S K 3  -  3  FNSKIft 
DO  417  I  -  HE0T1 .NELE.NSK3 
J1  '  I 

J2  I  +  2*NSK I  ft 

IFC J2.GT.NELE) J2  •  MELE 

URITEC6.222) CJ.UCJ) ,VCJ) ,YC(J) ,J -J1.J2, NSKIft) 

22  2  FORMAT  C 1 X , 3 < 1 5 , 1 X . 3E 1 2 . 5 , 2X >  ) 

417  CONTINUE 

1 F  C I FC . EG  *  0 ) GO  TO  520 
URITEC6.il) 

31  FORMAT (////. 42X.  49HVEL0CITY  FOR  SFLECTFD  VFLOCITY  STATI 
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IONS (ELEMENTS) ,111) 

10  010  I  1 , ISF  TE  ,  3 

1 1  —  I  +  1 

III  ■  1+2 

J  NSF'E(I) 

IF(  II  .GT  .  ISF'TEIGO  TO  016 
J1  '  NSFE(II) 

IF ( 1 1  I  .GT . ISPTE)GO  TO  017 
JO  =  NSPEdll) 

UNITE ( 6.222) J, U( J) r VC J)  ,YC(J> iJl >U(J1 > ,V(J1 > ,YC(J1 > ,J2,U<J2> ,V(J2) 

1 . YC  <  J2 ) 

00  TO  015 

516  WRITE  (6,223)  J.U(  J)  .V(  J)  i  YC(  J) 

GO  TO  015 

517  URITEC6.221 ) J , U ( J  > ,  V  (  J  ) ,YC(J) iJl.U(Jl) ,  V  (  J 1  > ,  Y  C  (  J 1  ) 

515  CONTINUE 

221  F  QRHAT (IX, 2(10, IX, 3E 12.0, 2X) ) 

22  3  FOR  MOT  (  1 X  ,  1  (  1 0  ,  IX  ,  3  E 1  2 . 0 , 2  X  )  ) 

2071  FORMAT  (9H  SOLUTION) 

0022  FORMAT  <4<1X,I4»1X.2F7.3>) 

2023  FORMAT (0< 2X, 1 1 , IX, F8 , 3)  ) 

23  FORMAT (/, 60X , 13HS0URCE  POINTS,/) 

24  F0RMAT( 1 X . 1 0 > OE 1 2 . 5 ) 

25  FORMAT(  1):,2<  I0.0E12.5)  ) 

520  IF  < ( IPC.NE . 0) .OR. ( I  PH AX . NE .0) ) CALL  CRUNCH ( NODFT r JFHAX , AF , U r V , NNQD » 
1MELE.HE0T  >XK.  IiONSP  ,  ISP  T  ,  I  PER  ,  I D  >  IPC  ,  HSFE  ,  ISPTE ,  IKT  ) 

99  CONTINUE 
STOP 
END 

*DECN  CF  UN 

SUBROUTINE  CRUNCH ( NODPT , IF'HAX , af , u, y , NNOD i NELE i nfot , XK  r IK , nsp  » ispt 
1 ,  I  PER,  III,  IPC.MSPE,  ISPTE  ,  IKT  ) 

c* *********************** ******************************************* ************ 

C  THIS  SUBROUTINE  PUNCHES  CARDS  FOR  AMPLIFICATION  AND  VELOCITY 

C* ********************************************************  ********************** 
DIMENSION  AF(NEOT) , U < NELE ) , V ( NELE ) , XK < 25 ) , I D ( 20 ) > MSP ( 1 00 ) , I  PER ( 5 ) » 

1AFT  (5)  ,  NSF'E  (  100)  >UT1  (3)  ,VT1  (3) 

URITE(6.101) 

101  FORMAT (  4  H  '  'C  ) 

WRITE <6, 10) ID 
10  FORMAT ( 20A4 ) 

IF  (IF’C.EQ.O)GO  TO  999 
URITE(6,50)ISPT, ISPTE, IKT 
50  FORMAT ( IX. 315) 

ISP  -  ISPT/5 
FISPT  FLOAT  (  ISFT) 

FISP  -  FISPT/5.0 

IF (FISP .OT . ISP  )  ISP  -  ISF  +  1 

I  COUNT  -  0 

J1  -  0 

J2  -  0 

HO  110  1= 1 , ISP 
ICOUNT  ICOUNT  +1 
J1  '  J2  +  1 
JO  Jl  +  4 

IF ( J2.GT . ISPT) J2  =  ISFT 
Ml  ■  0 

DO  112  J  -  J1.J2 
N 1  HI  +  1 
I.  =  NSP(J) 

112  AFT  (  PI  )  r  AF  <N  ) 

WRI TE ( 6. 11 1 ) ICOUNT , XK ( I K ) , ( AFT ( N ) , N - 1 . N 1 ) 

111  FORMAT  < IX , 15.F10 ,5.3F 13 ,4 ) 

110  CONTINUE 

ISPE  ISPTE/3 
FISF’TEsFLOAT  <  ISFTE) 
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FISPE=FISPTE/3.0 
IF  CFISFE.GT.  ISP  E>  ISF'E  =  ISFE  +  1 
K 1  0 

K2  0 
ICT  =  0 

PD  120  1=1  r  ISF'E 
ICT  ICT  +  1 
N1  -  K 2  t  1 
K2  K1  +  2 

IK  (  K 2  .  GT  >  I SF'TE  ) K2  ^  ISPTE 
Ml  0 

PCI  1  13  K- K1.K2 
Ml  =  Ml  +  1 
J  --  NSPECK) 

UTl(Hl)  =  -  U(J) 

113  U  T  1  (  M  1  )  ==  U  C  J  ) 

UNITE C 6.  11*.  IICT.XKCIK)  »  CUT1  CM)  »VT1  CM)  »K.  =1  .Ml ) 

11-1  FORNATC1X»I4>F10.3»6F11.4) 

120  CONTINUE 

UNITE ( 6 • ISO ) 

150  FORMAT (IX) 

WRITE  C  4  > 102  > 

102  FORMAT (  AH  ' 'L  ) 

I F ( I PM AX . NE . 0 ) GO  TO  1000 
RETURN 

99?  I F ( I PM AX . £0 . 0 ) GO  TO  2000 
1000  IC  0 

URITE(6.1010)C IFERC I > » Is 1.5) .XKCIK) 

1010  FORMAT C1X.5AA.F10. 5) 

PO  1100  1=1 . NNOP  r 11 
IC-IC  +  1 
LI  I 

L2"  L 1  +  10 

IF  <L2.GT.NN0D)L2  =■'  NNOD 
URITE(6i20)IC> (AF(L) . L=  LI . L2 > 

20  FORMAT! IX >IA .2X.11F4, 3) 

1100  CONTINUE 

PO  1200  I  =  1.NELE.5 
IC  -  IC  +  1 
LI  -  1 
L2  =  L 1  +  A 

IFCL2.GT.NELE)  L2  =  NELE 
URITEC6.30)  IC.  CUCL)  .  VCD  .L'=L1  >L2) 

30  FORMAT (IX. M . 2X  » 5 C 2F7 . 3  > ) 

1200  CONTINUE 

URITEC6. 150) 

2000  URITEC6.102) 

RETURN 

ENP 

»PECK  INF'T 

SUBROUTINE  INPUT2C  X.Y.NOD.NCON 
1 . NELE. NNOP. NOPR. NCS. NECT .NBANP.NCS2.PP. IPPC. PM. IPCH) 
c  «***M:***«*ftr ******************(******** **************************** 
C  THIS  REAPS  THE  INPUT  PATA  FOR  3-MOPF.S  TRIANGLE  ELEMENTS 

C**»  ******  *****  *************************  ******************  ************ 
COKiMON/PP/  UK  .  ALPHA  .  PI  .  PA2  .  PAT  .  RR  ,  UR  .  XL  .  WL  ,  RKHF  .  HKL 
PI  MENS  I  ON  XCNNOP) . Y < NNOP > . NOP ( NNOP > . NCON C NEL E . 3 ) t PP  <  NELE ) 
c* ********** *********************************************************** 
C  MOPC I)- EXTERNAL  NOPE  NUMBER.  IE.  ACTUAL  MOPE  NUMBER  ASSIGHEB 

C  X !  I  >  .  Y(I)--  X-ANP  Y-COORPINATES  OF  NOPCI) 

C  NCONCM.N)  li-TH-ELEKENT  NOPAL  CONNECT  I  V I T I ESE  N-I-J-K 

C*t ********>!  ********************************************************** 
N0PB1  -NN0P-N0PRF1 
N0PB2  NNOP-NOPR 
PO  18  I =N0PB1 >  NNOP 
IF ( I PP C . NE . 0 ) PPC I >=  PM 
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NOD(I) =1 

II  I -i;2Dbl  + 1 
N(V  II 

r.AD-  ( :  1-1 ) «  DA2 
y  (  I  >  -  RTOCOSCAAD) 

Y'I»-  RR4SIK'CAAD> 

is  continue 

CALL  CDEFTHC  I  DFC  i  DD ,  NNOD  »  N0DB2  •  NELE  ,  IDCH) 

IiD  310  I-  1.N0DB2.2 
I > -  1  +  1 

ir  C  I  .  EC.NODB2)  GO  TO  3314 

LEAD ( 4 .316)  N  0  D  (  1 > . X < I ) r Y ( I ) iNOIM  II).X(II) iY(II) 

GO  TO  315 

3314  READ<4»3316>  NOD  ( I  )  ,  X  ( I  )  »  V  (  I  ) 

313  CONTINUE 

IFUDCH.EO.  1) PRINT  31 6 , < NOD C I > . X < I > » Y < I >  1 1  -1  iNNOD) 

316  FORMAT (  2 <  1 1 0 . 2F 1 0 . 2  )  ) 

3316  FORMAT (  I10.2F10.2  ) 

I F ( IDFX  . ME . 0 ) GO  TO  400 
imii)  -  unci >/dm 

HO  60  I -"1  .NELE 
60  CONTINUE 
GO  TO  70 

400  I'O  56  1-  l.NNOD 

nnci)  -  ddid/dh 

56  CONTINUE 

70  no  415  L  =  1 iNELE.4 
LI-  L 
L2.  LI  +3 

1FCL2.GT.NELE)  L2  =  NELE 

415  F:EA  0(4.416)  (  (NCON  (LL.J)  » J-  1.3)  .LL-L1.L2) 

416  FORMAT (  4C5X.3I5)  ) 

IF ( IDCH. EQ, 1 ) PRINT  416. <  CNCONCL. J) » J =1 .3) .L =1 1 NELE) 

RETURN 

END 

FDECK  DEPTH 

SUBROUTINE  GDEPTH C I DFC . DD . NMOD . N00B2 : NELE . I DCH ) 

C***»tt*  **»*4:»<:«<!t»4(*»:**»****»  ***** ******* 

C  THIS  SUBROUTINE  READS  DEPTH  AT  EACH  GRID  MODE  FROM  CARDS  OR  USF.S 

C  SPECIAL  PROCEDURES  TO  READ  AVER.  DEPTH  FOR  EACH  ELEMENT 

C******  **********************  ***** ******* 

DIMENSION  DD ( NELE ) 

IF  < 1 DPC • ME • 0 ) GO  TO  1000 

C  SPECIAL  PROCEDURE  FOR  INPUT  OF  AVER  DEPTH  FOR  EACH  ELEMENT 

DO  100  I  -  1,121 
D D ( I >  -  24. 

100  CONTINUE 

DO  101  1=  122, 125 
DD ( I )  24. 

101  CONTINUE 

DO  102  I  126,318 
DD ( I )  =  24. 

102  CONTINUE 

DO  103  1.319.336 
DDCI)  =  24. 

103  CONTINUE 

HO  104  1=337.360 
DDCI)  =  29.5 

104  CONTINUE 

DO  105  I  361.371 
DDCI)  =  24. 

105  CONTINUE 

DO  106  I  372,385 
DDCI)  =  40. 

■106  CONTINUE 

DO  107  1-386,529 
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DD(I)=21 . 

107  CONTINUE 

HO  108  1-530 >2232 
JUKI)  -  38. 

108  CONTINUE 

no  109  I  -  2233 i 2333 
BD<I>  =■  42 . 

109  CONTINUE 
DIK2334)  =  38. 

IF(  IDCH.EQ.  1  )F'RINT  2010>  <  DC  <  I )  1 1  l.NELE) 

2010  FORMATt//. < 12F10.3) > 

GO  TO  2000 

C  READ  DEPTH  AT  EACH  NODE  FROM  CARDS  OF  PERM  FILE 

527  FORMAT  <  8 ( 15 iF5 . 0 ) ) 

1000  READ ( 4 . 527 ) < I ,DD< I > . X =1.N0DD2) 

IFdDCH.EO.l  ) PRINT  201  Or  <  DD  ( I  )  »  1--1  rNNOD) 

2000  RETURN 
END 

»0£CK  BAND 

SUBROUTINE  BAND (  NCON 
lrNELErNNODiNODRr NCSrNEQTtNBANDrNCS2> 

DIMENSION  NCON ( NELE  r  3 ) 

C**  *******  ********  ******  *********************  ************************* 

C  THIS  SUBROUTINE  DETERMINES  THE  BANDWIDTH  FOR  TRIANGLE  ELEMENTS 

C************************************************* *********** ********* 
MPS=  NODR+NCS 

kma:<=i 

DO  74  L!1 rNELE 
IKAX-NCON(L.l) 

IhIN- NCON(L> 1 ) 

DO  72  K-1.3 

IF (NCON ( L  r  K )  .LE.  IHIN)  I  HI N= NCON < L . K ) 

IF ( NCON <  L » K)  .GE.  INAX)  IMAX-  NCON(LrK) 

72  CONTINUE 

KK-IMAX-1MIN 
IF (KK  .GT.  KM AX)  KHAX-KK 
74  CONTINUE 

NBAND-KMAX+1 

NBAND-  AHAXO  <  MB  AND  >  MF'S ) 

WRITE < 6.84 )NB AND »HB AND »MPS 

84  FORMAT  < // //1X»40(1 HO ) » /  I2H  BAND  WIDTH*  r 13. /27H  BAND  WIDTH  W.R.T.  E 
lLEMENT*rI0r/34H  BAND  WIDTH  W.R.T.  SOURCE  +  BDARY-rI5i/lXr40<lH0)r/ 
2) 

RETURN 

END 

*DECK  ASBK 

SUBROUTINE  ASF MBK ( X r Y r NCON t SVSK r SVSQ 
1  rSTrMGL 

If NELE .NNOD. NODR . NCS.NECTr NBAND >NCS2rDD i  NELE3 

1  . NLGG . I DPC  r I  BUG ) 

C* ******** ************************ ************************************ 

C  THIS  SUBROUTINE  ASSEMBLES  ELK ( I r  J) - ELEMENT - INTO  SYSK(I.J)  — 

C  -  SYSTEM 

C  ELK-  ELEMENT  STIFFNESS  MATRIX  3  BY  3 

C  SYSK-  SYSTEM  STIFFNESS  MATRIX  NEOT  BY  NBAND 

C*****»( ****************************************************** ******** 
COMPLEX  SYSKrSYSOrST 

COMMON/DB/WKr ALPHA. D1 r DA2 . DA4 r RR . WP . XL r WL r RKHP . HKL 
DIMENSION  X(NNOD) . Y < NNOD ). NCON < NELE . 3 ) .SYSO(NEOT) 

COMMON  /NEW/  NUMBLK.MLGiNL0L.Lll.LQ.L4.L3 
DIMENSION  NR <  3 ) . ELK < 3 . 3 ) . DD < NELE ) 

DIMENSION  NR ( 3 ) 

DIMENSION  ST(NGL. NBAND) 

C  CHANGES 

DIMENSION  SYSKCNLGG. NBAND) 

LEVEL  2  >  SYSK  >  SYSO . ST 
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DO  122  I-l.NEQT 

122  EVSQ (I)--(O.OtO.O) 

L2  0 

C  CHANGES 

L3-  NLGiNBAKD*2 
LO=MLG 12 

do  eee  iq- i.numblk 

LEM  NLG 

IF  (10  .EQ.  NUHBLK )  LEN-- NLGL 
L2-L2+LEN 
LI  (  IC-1  MNLG+1 
LU-  Ll-1 
I'O  12  I  I.  LEN 
C  CHANGES  MADE 

DO  12  J  l.NBAND 

IF  <10  .ME.  NUMBLK)  GO  TO  7176 

ST(I.J)  -  ( 0 . 0  >  0 . 0 ) 

GO  TO  12 
7176  CONTINUE 

SYSK(I.J).  <0.0. 0.0) 

12  CONTINUE 

DO  25  L-l.NELE 
IGO  =1 

DO  14  J= 1 , 3 
MR<  J) =0 

NR ( J ) =  NCON  <  L  >  J ) 

IF  < (NR< J) .LT.L1 ) .OR. <NR< J> .GT.L2) >  GO  TO  14 
IGO-2 
MR  <  J  )  -  1 
14  CONTINUE 

IF(IDUG.EO.O)GO  TO  1000 
URITE<6.10)(NC0N(L,J),J-1>3) 

10  F0RMAT(1X,4H**»*.3I10.4H*»**) 

1000  GO  TO  (25.44) »IGO 
44  I1-*NR(1) 


I2  =  NR  <  2 ) 

I3-NR ( 3 ) 

XDX(Il) 

X2-X<12) 

X3=  X< 13) 

Y1 =Y< II) 

Y2-Y  < 12 ) 

Y3-YCI3) 

1 F  < I DPC • NE . 0 )  GO  TO  210 
D?  =  D D  <  L ) 

GO  TO  200 
210  DID DD(Il) 

D 1 2 - D  D  < 12) 

D13-  DIKI3) 


D9=(Dll+D12+D13)/3. 

IF  <L.GE.  1435.AND.L.LE. 1571) D 9=0.5 

C*»wmtmt  SPECIAL  STATEMENTS  HAY  BE  ADDED  TO  COMPENSATE  FOR  ****»»«** 

mmmtm 
mmm 
NELE3 ) 

i  r  v  i  puu  .  cu  .  <J  ;  uu  iu  iuju 

WRITE <  6. 15) < (ELK(H.N) .N  -K. 3) «K  -1 . 3) 

15  FORMAT (6E15.6) 

1050  DO  20  1-1.3 
DO  18  J-I.3 

I F  <  NR (J)-NR(I)  .GE.  0  )  GO  TO  16 
IF  ( HR ( J )  .EC.  0)  GO  TO  18 


C************  AVERAGING  PROCEDURE.  FOR  EXAMPLE.  AS  ABOVE-- 
C*»*  IF<L.GE.  1435.  AND.  L.LE.  1571)  D9  0.5 
C*»**»:***»»»:***»** 

200  DIN  -•  1.0 

CALL  ELMKIXl .X2.X3. Y1 , Y2.Y3.WK. D9i DIN, El.K .AREA »L» 
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LR=NR(I)-NR(J)+1 

J1  =I\'R  (  J )  -L 1 1 

IF  <10  .NE.  NUNBLK)  GO  TO  7177 
ST(J1,LR)=ST<J1*LR)+ELK(I,J) 

GO  TO  18 

7177  CONTINUE 

SYSK(JltLR) =  SYSK ( J1 > LR  HELK< I » J) 

GO  TO  18 

16  IF  ( MR ( I )  . EQ.  0)  GO  TO  18 
LS  NR<  J)-NRU>  +  1 
J2  =  NR ( I ) -LI 1 

IF  (IQ  .NE.  NUNBLK)  GO  TO  7178 
ST(J2,LS)  -ST  (  J2  > LS) f ELK < I » J ) 

GO  TO  18 

7178  CONTINUE 

SYSK< J2,LS>=  S)SK< J2,LS>  +  ELK(I» J) 

19  CONTINUE 

20  CONTINUE 
25  CONTINUE 

IF  <10  , EQ .  NUNBLK )  GO  TO  888 
L4-I0+10 

CALL  ECUR(SYSKtL4rL3> IERR) 

IF ( IERR.EO. 1 )G0  TO  999 
888  CONTINUE 
54T  RETURN 
999  PRINT  777 

777  F ORM  AT  <  20H  WRITE  ERROR  IN  ABSK  ) 

STOP 

END 

♦DECK  GSPK 

SUBROUTINE  GSPK < SYSK t X J, XV r XH . DH. ANGL 
1 >NELE»NN0D,N0PR,NCS«N£GT»NBAND>NCS2,N0PR1, VJ,VY»NLCG,IBUG) 

C  CHANGES  HADE 

C*t*l****t****»***l»t*t*t«*t»tt*tttl4»************(*S*S*******«»t**«*t 
C  THIS  SUBROUTINE  OBTAINS  MATRIX  K2  +  K3,  AND  STORES  IN  THEIR 

C  CORRESPONDING  ADDRESSES  IN  SYSKU.J) 

C********»:»:M:»»**»**t  *«»:»♦********»******»***»****»**♦*************»* 

COMPLEX  XH i DH i TH  j  SYSK 

COMMON/DB/  UK, ALPHA, D1 ,  PA2  ,  DA4  ,  RR ,  UR ,  XL  »  UL  .  RKHP .  HKL 
DIMENSION  SYSK  (NLGG ,  NB AND  )  ,  X J(NCS) , XY <NCS >  , XH( NCS  >  >  DH(NCS) 

1 , ANGL ( NODR  1 ) 

DIMENSION  VJ<1),VY<1> 

LEVEL  2, SYSK 

COMMON  /NEW/  NUMBLK,NLG,NLGL»L11 ,LS 
1 ,LA,L3 

C*»:*T*»*  ******»»'*»»:***  I-***************:**********************  **«»♦****« 

C  OBTAIN  BESSEL  FCN  XJ.XY.XHiDH 

C  *«**«»*  *4:  ****:*:*:»  t4:**t*«f***«:*|:»*»4:*»t  *******  t  ***««*«tt*«******t  ****** 
UR-  UK*RR 
XL-  RR4DA2 
UL  UKTXL 
do  12  r  1 , NCS 
11=1-1 
FNU  0. 

CALL  BESSY<UR, FNU, II , VJ, VY ) 

XJ< I)  V J ( I ) 

X  Y  < I >  VY ( I  ) 

I F ( I  BUG . EQ • 0 ) GO  TO  12 
WRITE<6»20)  XJ(I),XY<I) 

12  CONTINUE 

DO  15  I  -I, NCS 

XN  =  FLOAT< 1-1 )/UR 

XH< I  )-  XJ< I >  +  <0.0, 1 .0)*XY(I ) 

IF  <  I  .NE.  1 )  GO  TO  9 

DH< 1 )  -XJ<2>-< 0.0,1 .0)»XY<2> 

GO  TO  10 
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9  DH(I)*(XJ(I-1)-XN*XJ(I) )+(Q . 0, 1 . 0) *(XY(I-1)-XN*XY (I) ) 

10  CONTINUE 

IF<  IBUO.EQ.O)GO  TO  IS 
WRI T£ ( 6  *20 )  :<H,I)*DH(I) 

15  CONTINUE 

20  FORMAT (OH  OSPK  BY- » 4E15. 6) 

1520  FORMAT  (20H  BLOCK  ERROR  IN  GSFK  .15) 

C  STORE  THE  UPPER  TRIANGLE  ELEMENT  OF  K2  IN  CORRESPONDING  ADDRESSES 
C  IN  SYSK(I.J) 

cimutmmmummmmttmmtmmmmmmmmmm 

RKHP-  3.14159265*UR*D1 
DO  38  I-l.NCS 
II ’NN0D+I-L11 
NLLL=1 

IF  (Il.LE.O  .OR.  Il.GT.NLGG)  GO  TO  15S1 
CX- 0.5 

IF(I.EG.l)  CX  1.0 

StSK<Il.l)--SYSK(Il.l)+CX*RKHP«:XH(I  ,*DH(I  ) 

38  CONTINUE 

cmmmmmummmmummmmmmmmmmmtm* 

C  STORE  K3  IN  CORRESPONDING  ADDRESS  IN  SYSK(I.J) 

Dmmmmmmmmmmtmutmmmtttmtmmmmtm 

DO  40  I  = 1 . NOD R 1 
40  ANGL(I)- < 1-1 )*DA2+DA4 
HKL  -WLTD1/2 . 

N1-NN0D-N0DR 
DO  48  I  -1* NODR 
II-N1+I-L11 
NLLL=2 

IF  (II.LE.O  .OR.  Il.GT.NLGG)  GO  TO  1551 
n-i-i 

NODR  I -NO DR +  1-1 
DO  48  J’l.NCS 
JJ=N0DRI+J 

IF <  J  .GT.  1)  GO  TO  42 
TM  HKL4DH ( 1 ) 

IF  (I  .GT.l .  AND.  1  .LT  .NODR)  T  H  ■  2 . 0  ♦  T  M 
GO  TO  46  , 

42  J2--J-1 

IF(I.LE.l.OR.I.GE.NODR)  GO  TO  51 
A1  -ANGL ( I ) 4 J2 
A2«ANGL< II )»J2 

TM*HKL<DH<J)t(C0S(A2)+C0S(Al>) 

GO  TO  46 

51  IF(I.EQ.NODR)  GO  TO  52 
A1  ANGL(I)*J2 
TM*HKL4DH< J)»C0S(A1> 

GO  TO  46 

52  A2;  ANGL  < 1 1 ) * J2 
TM=HKL»DH< J)*C0S<A2> 

46  SYSKdl,  JJ).  SYSK(II.JJ)-TM 

43  CONTINUE 

IF  <  NUMBLK  .EQ.  1)  RETURN 
L5  'MLGG4N8AN042 
L6*L  4  + 1 

CALL  ECUR(SYSK.L6.L5. IERR) 

IF(  IERR.EC. 1 ) GO  TO  2000 
IF(  IBUO.EO.OGO  TO  3000 

PRINT  1018.  <SYSK(I.1).SYSK<1.2).I:  1.NLGG.4) 

1018  FORMAT  (IX. 4014. 3) 

3000  RETURN 
2000  PRINT  999 

999  FORMAT  (20H  WRITE  ERROR  IN  GSFK  ) 

STOP 

1551  PRINT  1520.  NLLL 
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* t ECK  LOAD 

SUBROUTINE  LDAD<SHfSYSCfXJfDHf ANGL 
IfNELEfNNODfNODRfNCSfNEQT  fNBAMDfNCS2iNDDR1 f I  BUG) 

C  THIS  SUBROUTINE  GIVES  LOADING  TERNS  6  SYSQ(I) 

C*«******X**«**<*«  **lt******t***************************************** 
COMPLEX  DUMY . SM . SYSQ  f DH . TN  >  CM . DUH2 

COMNON/BB/  UK  >  ALPHA  r  D1 1 DA2 f  DA4 » RR  f  UR  r  XL  f  WL  » RKHP  fHKL 
DIMENSION  SM <  NQDR1 ) >  SYSC ( NEOT )  f  DH  < NCS ) f  X  J  ( NCS ) t ANGL ( NODR1 > 

LEVEL  2 f SYSQ 
DULu 0. StDIXUL 

C********************************************************************* 

C  CALCULATE  LOADING  TERMS  U.R.T.  NODAL  POTENTIAL  +  STORE  IN  SK(I> 

C  *<:X*T.t*X*  *****X*|:  **************************************** ************ 
DO  12  I-lfHOORl 
THE  T  A-  ANGL ( I ) 

CX-COS< THETA- ALPHA) 

AUG-WRXCX 

DUMY=  COS( AUG)  +  (O.Of 1 .0)*SIN<AUG> 

C2=C0S<TH£TA+ ALPHA) 

AU2-WRXC2 

DUM2-COS<AU2)+<O.Oi 1 .0)*SIN(AU2> 

SM  ( I )  --  DUL*  ( 0 . 0  f  1 . 0)  *  <CX*0UKY+C2*BUM2> 

12  CONTINUE 

IF ( I  BUG  .  EC . 0 ) GO  TO  100 
WRITE (6f 3) <SM<I) fI-IfNODR) 

8  F0RHAT(//f<7E15.7)> 

100  N 1  -  MNOD-NODR 
DO  16  I‘*  If  NODR 
I1=N1+I 

IFU.EQ.l)  GO  TO  1 A 
IF< I  .EO.NODR)  GO  TO  IS 
1 1  •  I  - 1 

SYSQ(I1)-SYSQ<I1)  +SK< I ) +SH( II) 

GO  TO  16 

14  SYSQ<I1)*SYSQ(I1)+SK(1> 

GO  TO  16 

15  SYSQ(I1)*SYS0(I1)+SM(N0DR1) 

16  CONTINUE 

IF  < I  BUG . EQ ■ 0 ) GO  TO  110 

URITE(6f8> (SYSC(I>fI~N1fNN0D>  . 

C  »«******** ********* ****** **************************************** 
C  CALCULATE  LOADING  TERMS  U.R.T.  SOURCES  AND  STORE  IN  SYSQ(I) 

c*mt*m*tmm**t*m««*«**»«*****»*****”**”**m**MW*MM* 
110  RKHP2  -2.0XRKHP 
DO  40  I - 1 f NCS 
NNI=  NNOD  +  I 
II--I-1 
A I -ALPHA* I 1 
CM;  <0.0# 1.0>**II 
T  M  *  C  M  $  X  J  < I ) t  DH ( I >  *  C  OS ( A I ) 

SYSOINNI )‘  SYSO(NNI )-RKHP2«TM 
40  CONTINUE 
RETURN 
END 

* D E C K  ELMT 

SUBROUTINE  ELMK ( XI f  X2 1 X3  f  Y1 f  Y2f Y3fUKf  Df  D I N f  ELK f  ARE A f  LUM f 

C4t **»  »*»»»**»*« »***»*»« ***  *********  ********************************** 
C  GENERATION  OF  TRIANGLE  ELEMENT  MATRIX  ELK 

C  BIN-1.0  IF  IN  OUTER  DOMAIN?  D'DIf  HELMOTZ  EO. 

C  DIN=0.0  IF  UNDER  FLOATING  BODY1'  D=D2 f  LAPLACE  EQ. 

C» ***<***»» *********************************************************** 
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DIMENSION  ELK(3,3) 

DO  8  1=1.3 
DO  8  J'  1 » 3 
8  ELK ( I »  J  > -  0. 

81-  Y3-Y2 
B2-  Y1-Y3 
B3-  Y2-Y1 
Cl*  X2-X3 
C2-  X3-X1 
C3-  >;i-X2 

i*«REft  0.5t(81*C2-B2*21) 

G  =  32 . 2 

FAC  -2.*AREA/G 
LUH 1 - 3*LUM-2 
LUH2-3*LUM-1 
LUM3=  3»LUM 

IF <  AREA  .Gf.  0.  )  GO  TO  10 
URITE<6»100>  AREA.LUM 

100  FORMAT ( IX i 23H  DEBUG  ELMK  100.  AREA  =.E12.5»5H  AT.I5.11H-TH  ELEME 
1NT ) 

STOP 

10  A4-4.*AREA 
UKVD-UK 
UN~UK/SQRT  <  D ) 

WKA=  DINTWKtUKT.  AREA/12. 

ELK ( 1.1)=  (  (Bl*Bl  +  CltCl)/ft'»-2.TWKA>»D 
ELK ( 1.2)*  ((BHB2+C,.  C2J/A4-  UKA>*D 
ELK  < 1 . 3 ) =  < <B1»B3+Ci <C3)/A4-  WKA)*D 
C  ELK(2.1)-ELK<1.2) 

ELK (2.2)-  <  (B2*B2+C2J'C2>/A4-2.*UKA>»D 
ELK (2.3)"  ( <B2*B3+C2»C3)/A4-  UKA)*D 
ELK<3.1)-£LKU»3> 

ELK(3.2>-ELK(2.3) 

ELK ( 3  »  3 ) *  ( <B3*B3+C3*C3>/A4-2.»UKA>»D 
UK -UK  VD 
RETURN 
END 
•DECK  BEST 

SUBROUTINE  BESSY (X.FNU.N.BJ.BY) 

C  BESSY  NYU  MATH  UTILITY  SUBROUTINE  FOR  BESSEL  Y  3/13/64  F.  RAGUSA 

DIMENSION  BJ(l).BYCl) 

C  DIMENSION  B J ( 1 ) . BY ( 1 > 

EQUIVALENCE (FN.NF) 

C 

XSAVE =X 
M1  =  N 

MF- 131070 
B  J  ( 1 )  s:Fh 

CALL  BESSJ(X.FNU.FN.BJ) 

BY  ( 1 )  1:FN 
YNU- BY ( 1 ) 

N-Hl 
BY  <  2 ) *X 
X*XS AVE 
C 

NN= IABS(N) 

N1--MN-1 
CONST -2 • 0/X 
F  I-.31415926SE  +  01 
C 

IF(»10.0)  30.19.19 
30  X 1 0  =X  +  25 . 0 

kio  - >r i o 

N10-NN+10 
M»  HAXO ( K10.N10) 

C 
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IF  (X  -  1.0)  1,  1,  3 

1  KF'-  172.69386/ (3. 6888795- AL0G(X>  > 

M  =  i1IN0<M,KF> 

GO  TO  2 

3  KP= 39. 0*X4*. 3333333 
M  =  rlIN0<M,KP> 

2  M  li/2 
K=2*M+1 
ARG  FNU*PI 

GARG=-GAMMA<  1 .  0  +  FNU.0.  >**2 

COMPUTE  GAMMA  ZERO  ECU.  15, PG.  0  t  IP  KU-0  USE  EQU.  16. 

16  IFCFNU)  10.16,15 

15  TERM  M1.0/PI)»  CONST  **<2.0*FNU> 

GAH1  -  COS ( ARG  > /SIN ( ARG >  -  TERM* (GARG/FNU) 

COMPUTE  GAMMA  ONE  ^OR  K‘U  NOT  EQUAL  ZERO. 

GAN2- 2 . 0  *  TERM  *  GARG  *  (FMU+2.0)  /  (l.O-FHU) 

GO  TO  10 

COMPUTE  GAMMA  ZERO  AND  ONE  FOR  NU  EQUAL  ZERO. 

16  TLOG  •-  ALOG ( X/2.0) 

A=.577215664E+00 
PIH =2. O/PI 

GAM1=  PIH< ( A+TLOG) 

GAM2- 6 . 0/P I 

COMPUTE  GAMMA,  YNU,  AND  BY<1>»  AND  BY  (2) 

EQUATIONS  15,17,  AMD  18. 

10  BY (2)  =0.0 
BY ( 1)= 0.0 
GAM3 =0.0 

E1--<1.0/PI)|:C0NST  **< 1 .0+2.0*FNU>  *  GARG 
BY<2>  E1*BJ<1>+BY<2> 

El=  GAM1-GAH2/2.0 
B Y ( 2 )  :E  1  *B  J  <  2 )  +  B Y  <  2 ) 

YNU=GAH1»BJ<1> 

TXNU=  3.0  +  FNU/X 
AB- ABS  <  B J ( 1 ) )  -  .000005 
12  =  1 
MP1-H+1 


DO  11 
12  =  I 
FI  =  I 
FIM  = 

F  1 2  - 
DENOH 
GAM3  = 
GAM3  -= 
YNU  = 
IF< AB) 
WHEN  J 


I-2.MF1 
2  +  2 


=-  FI  *  <  F I  -FNU  >  *  <  FNU+FI 2-2 . 0 ) 

( FNU+F 12 ) t ( 2 • OtFHU+F IM ) * ( FNU+F IK  > /DENOH 
-GAM3*GAM2 
GAM2*BJ< I2)+YNU 
18,18,23 

(NU)  IS  NEAR  ZERO  COMPUTE  BY<2>  FROM  EC.  18, ELSE  EQ.  17. 


E1-TXNU1GAM2 

BY<2>  E1*BJ< 12  >  +  BY  C  2 ) 

IF ( 12-K )  25.110,110 
El  --<GAM2-GAM3>/2.0 
B Y ( 2  )  =E1*BJ< 12+1 >+BY(2> 

GAN  1  =  GMM2 
3AM2= GAM3 

BY ( 1 ) =  YNU 
I  +  <  A  8 )  19,19,17 

BY(2)  =  (  YNU+BJ(2)-2.0/(PI*>:>  )/BJ(  1  > 
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o  n  n 


C  IF  N-0  OR  1  GO  OUT  ALL  BY  S  CONFUTED . 

19  I F ( N )  30.49.20 

20  IF(N-l)  49.49.21 

CONFUTE  Y ( N+ 1 )  BY  RECURRENCE 

21  HO  22  I- 1 ,N1 
F 1 1  =  1 

22  BY< I+2)  =  CONSTt(FIl+FNU)*BY( I  +  1)-BY(I) 

49  RETURN 

50  IF (FNU-  .500000)  54.53.54 

53  BY  <  2 ) -B J ( 1 ) 

IF<NN-1)  49.61 .60 

60  ARC=-t.O 
B  Y  (  3  ) :  ARGtB J  <  2  ) 

BJ(2)-C0NSTtFNU*BJCl >-BJ<2> 

ARC =  -ARG 

IF  (NN-2)  49.51.56 
56  HO  55  1=2.  N1 

BY< l+2>- ARG*BJ< I+l) 

55  ARG=-ARG 
GO  TO  51 

54  BY<2)= CQNST*FNU*BYC1)-BY<2> 

BJ<2)=  CONST*FNU*BJ<  1  )-"BJ(2) 

IF(NN-t)  49.49.51 

51  FRAC=  FNU 
HO  52  I-1.N1 
FRAC=-  FRAC-1  •  0 

BJ<  I+2)«-C0NSTtFRACtBJ(I  +  l)-BJ(I> 
IF(FHU-. 50000)  58.52.58 
58  BY(I  +  2)=  CONSTtFRAC».BY(I  +  l)-BY(  I) 

52  CONTINUE 
GO  TO  49 

61  BJ(2)-C0NST*FHU<:BJ(1)-BJ(2) 

GO  TO  49 

ENH 
*HECK  GAHA 

FUNCTION  GANNA(U.U) 

DIMENSION  B ( 8 ) 

BC1)  =  .035868343 
B(2)=-. 193527818 
B  (  3  )  •=  .482199394 
K 4 >=-.756704078 
B<5>*=  .918206857 
B<6)  =-.897056937 
B  (  7 )  -  .988205891 
K  8 )  =  -  .577191652 
A  =  U 
X  =  U 

XN  =13.0 

1  IF  <A>4.2.3 

2  IFCX-1 .0)4000.4000.5000 

4  2  =  A I  NT ( A) 

Y  =  ABS<  Z  -  A  ) 

EF'S  =  3  .  OE  -  8 

IF  <Y-EFS>21»21»5 

21  A  =*  Z 

GO  TO  22 

5  Y  1 . 0-Y 

IF  < Y-EFS)6.6.3 

6  A-Z-l .0 

22  I F  <  X- 1 .0)4000.4000.5000 

3  IF(X)7. 1000.7 

7  IF ( X-l .0)9. 9. 8 
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8  IF  (A)  5000,  10,  10 

10  IF< (X/A)-l .0)1000*3000*5000 

9  IF  <A&S<A>  -1.)  11.1000.1000 

11  IF<X-.  1  )1000. 12*12 

12  JF(X-.2)13*14*14 

13  XN-140.0 
GO  TO  0000 

14  IF<X-.4)15*16*14 

15  XN  *  80 . 0 

GO  TO  5000 

16  IF ( X- .  6)  17 . 18, 18 

17  XN*  60  •  0 

GO  TO  5000 

18  IF  <  X- . 8  > 19  *  20*  20 

19  XN-40.0 

GO  TO  5000 

20  XN*  20.0 

5000  XFT *X 
NxXN 

DO  5001  I  *  1  *  N 
XFT- (XH/XFT) +1 .0 
XFT-< (XN-Al/XFTJ+X 

5001  XN-=XN-1.0 

TEH  x  <ALOG<X)»A)-X 
GO  TO  6000 

4000  SUM=0.0 
ANS  *XFT 

XFT  =  £XP<  TEH)/XFT 

EPS  =1  .OE-8 

XM-1 .0 

XhT  =  1,0 

EX**- 1  , o 

TEM  =  X 

4001  Y* TEM/(XNT*XN) 

4011  SUM-SUH+EX4Y 

IF  <  ASS ( Y ) -EPS )  4003  >4002*4002 

4002  TEMxTENFX 
Xh=  XH  +  1 .0 
XMT -XMT4XH 
EX* -EX 

GO  TO  4001 

4003  E* -ALQG<X)-.57721566-SUH 
IF  <A)4005r4004>  4005 

4004  ANS*E 

GO  TO  6000 

4005  IF<A+2. 0)4010. 4009*4008 

4008  ANS*-E+EXP<-X)/X 
GO  TO  6000 

4009  sun-<i.o/x)-<i.o/<x*x)> 
ANS-*<E-EXP(-X)*SUI1>».3 
GO  TO  6000 

4010  EX  *- 1 . 0 
N— A-l  .0 
XMT*  1 . 0 
SUH*0 . 0 

te;i*x 

DO  4006  1*1. N 
TEM  *TEM*X 
XM**I  +  1 
Y  =  XI'.T/TEM 
SUM -SUH+EX4Y 
EX  -EX 

4006  XMT*XMT*XM 
SUH-SUJI+l.O/X 
2*=E-EXP(-X)*SUM 
Y* A  8S ( A ) 
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EX— 1.0 
XflT-  X  .  0 
N  -  Y 

4007  1  =  2. N 
Xh=I 

XMT  =  Xf.TYXh 
4007  EX-  -EX 

ANS =ZLEX/XMT 
GO  TO  6000 

1000  IFIA-1. 0)1001. 1003.1004 

1004  IFIA-2 .0) 1003 . 1003, 1002 
1003  AT«A-1.0 

ATM  =B< 1 >*AT 
DO  1005  1=2.8 

1005  ATM=(ATK+B(I>>*AT 
AN  j-ATM  +  1 .0 

GO  TO  2000 

1001  IF<A>1007, 1006. 1006 

1006  BN--A 
BS  =  A 

GO  TO  1010 

1007  BN  =  A 
BS  =  A 

1011  IF(BS) 1009, 1010, 1010 

1009  BS=BS+l.O 
BN= BS*BN 
GO  TO  1011 

1010  ATM=BS*B(1> 

BO  1012  1=2,8 

1012  ATM  =  < ATh  +  B < I ) ) TBS 
ANS =<ATM+1 .0)/BN 
GO  TO  2000 

1002  BN=A-1.0 
BS=  BN 

1013  IF<2.0~BS)1015, 1014.1014 

1015  BS=BS-1.0 
BN=BS*BN 
GO  TO  1013 

1014  BS=6S-1.0 
ATM--BS*B<1) 

BO  1016  1=2.8 

1016  ATM=<ATK+B(I))*BS 
ANS =  (ATHU ,0)*BN 

2000  IF<X)2001 .6000,2001 

2001  W-1.0 
XM=  1 .0 
SUM=  1  .0 

E  F'  5  -  1  .OE-8 

2002  W  (X/(A+XN))*U 
SUM =SUH+U 

XM=  XM+1 .0 

IF  ( ABS(U)-EF'S)  2003,2002.2002 

2003  T=A*ALOG(X> 

T=EXF(T-X)/A 
ANS  -ANS-(T*SUM> 

6000  GAMHA= ANS 
RETURN 
END 

*  D  E  C  K  BESJ 

SUBROUTINE  BESS J < X , FNU , FN » F ) 

C  BESS J  NYU  MATH  UTILITY  BESSEL  J  SUBROUTINE  3/15/64  F 

DIMENSION  F ( 1 ) 

EQUIVALENCES, MF) 

C 

N  =  FN 

NN=IABS(N> 


.  RAGUSA 
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C  WHEN  BESJ  HAS  BEEN  CALLED  BY  BESY  F(l)  WILL  CONTAIN  INTERCER  131070 
FH=F(1 ) 

I Y1  1 

IF  ( KF- 131070 )  2.1.2 

1  IY==2 

2  CONST-2 . 0/X 
C 

IF(X-GO.O)  6 1 8  i  8 
6  :<10-X  +  23.0 

k  i  o  -  >:  i  o 

N 1 0  NII  +  10 

40  M-r.AXOCKlO.K'lO) 

C 

IFCX-l.O)  3 > 3 1  4 1 

3  K  P '-  (  172. 69388/(3. 6888795- A.  LOG(X)  )  )  -1  ,0 
N=MIN0(M.KP) 

GO  TO  4 

41  KP= 39. 0*X*t. 3333333 
M  =  M I NO ( N »  KP ) 

4  M  K/2 
43  K-24IU1 

K2=  KF1 
K3=  K  +  2 

F  <  K2  > • 1.0E-3S 
F  ( 1.3 )  •  0 . 0 
C 

HOC  L*1 »K 

I  =  R  +  1  -L 
F  L I  =  I 

5  F  (  I  >  -CONSTKFLI+FNU)  »  F(I+l)-F(I+2> 

C 

C  F I  Nil  ALPHA  EITHER  FROM  E2UA.  12. PAGE  4  UHEN  X  IS  LESS  THAN  8.0 

C  OR  FROM  VESSEL  FUNCT.  FOR  LARGE  ARGUMENTS  WHEN  X  IS  1 OR  GREATER. 

C 

IF(X-IO.O)  7.8,8 
7  PHI  - FNU+2 .0 

ALF  =F'HI*F(3)4F  (  1  ) 

M0=  3 
C 

noio  1=2.8 
MO- MO +2 
FM2  =2*1 
FMl  1-1 
F  I  =  I 

TEMF-( (FNU+FK2 )»  (FNUFFNI )  > / ( F I  * ( FNU4FU2- 2 • 0 ) )*FHI 
F HI -TEMP 

15  At.F1  PH  I  *F  (HO  )  +ALF 
C 

ALF  CONST *>  FNU* GAMMA ( FNU  + 1 .0.0. ) * ALF 
C 

II  =  1 
C 

C  F I  Nil  J  (  N  )  EQUATION  7, PAGE  3. WHEN  X  LESS  THAW  10.0 

C  OR  J  <  2  >»  J ( 3 )»...>  J ( N )  F ( 2  > /ALPHA ..... F ( K ) /Al PHA  WHERE 

C  ALPHA- F  <  1  ) /J<  1  )  AND  JU>  A*COS<FHl>  FROM  PATH  8. 

C 

16  0017  I  II rK 

17  F  < I >  =F  <  I  )  / ALF 

C  TEST  FOR  NEGATIVE  N  AND  RECOMPUTE  FS  BY  RECURSION. 

IF(H>  18.22.22 

18  IF ( I  Y  —  2 )  19.22,19 

19  F < 2 >  -CONST*FNU»F( 1 ) -F (2) 

IF(l.'N-l)  22,22.20 

20  FRAC- FMU 
N 1  K'N+1 

00  21  L-3.N1 
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FRAOFRAC-l.O 

F  i  L  )  -  COI!ST*FRAC*F  ( L  -  1 )  -  F  (L-2 ) 

IF  ENTRY  WAS  FROM  BFSY  RETURN  A  VALUE  FOR  Y<1>  IN 
RLACE  OF  N.IY  WILL  RE  ECUAL  TO  2  IN  THIS  CASE. 

22  GO  TO  (O'). 23)  .  IY 

2  3  FN  Y 1 
X  =  Y2 

2-1  RETURN 

RATH  6.  WHEN  X  IS  GREATER  THAN  OR  -TO  10.0 
CONFUTE  J(l)  FRON  II. GOLDSTEIN  PARER  BESSEL 
FUNCTIONS  FOR  LARGE  ARGUMENTS. 

8  FOUNT-  1 
GNU=  FNU 
CO  ■ .25 
Cl  =.15425 
C2  -  - . 375 
C 3  -.1171875 
C4  =-1.15425 
C5  1.875 
C6  =• .952148437E-01 
C7  -=-2.38671875 
C8  = 14 . 2265625 
C9  =---19.68750 
CIO  =  .809326171E-01 
Cl  1  -  4 1 0058593E+01 

C 1 2- . 582246093E+02 
C13  =-277.87500 
014=  354.3750 
C15 - . 416666666E-01 
Cl 6--  -  .25 
C 1  7 ==  •  125E-1 

C19==~558035718E-03 
C20--  . 424 1 07 1 42E+00 
C21=  3.60267857 
C22  =-5.625 
C23-.30381944E-02 
024  =-.486111 
C25=  .  102864583E  +  02 
C26- -58.0 
C27  =78.75 
9  AL 1 = GNU** 2-. 25 
A2  -=C0*AL1 
A  4  =C1*AL1 
A4-(A4  +  C2)>  AL 1 
A6=C3*AL1 
A6=(A6+C4) * AL 1 
A  6  = ( A6-FC5) *AL1 
A8-C64AL1 
A3=  <  A8  +  C7  )*ALl 
A8=(A8+C8  >*AL1 
A8- ( A8RC9  >  * AL 1 
A10  C10YAL1 
A10  -  < A10  +  C1 1 ) *AL1 
A10-( A10tC12)*ALl 
A10  ( A 1 OtC 1 3 ) » AL 1 

A10  ( A 1 0  +  C 1 4 ) * AL 1 

FI- .314159265FF01 
TS- 1 .0/X 
T2=-TS*«2 
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B=A10*T2+AB 
B=B*T2+A6 
B-B4T2+A4 
B=fi*T2  +  A2 
BNU  -BJT2+1 .0 

ANU  BNU/BCRT  (  .S»FI*X> 

RAGE  20 1  EQUATION  12  TO  GET  F'HI  ZERO 

A2  = . 0*  AL 1 
A4  =C1G#AL1 
A4  =<A4+C16)tALl 
A6  =C17*AL1 
A6-<  A64C18MAL1 
A6- ( A&  + . 7G)*AL1 
A8  ;C  1  9# AL 1 
A8=( A8+C20)*AL1 
A8=< A8+C21 ) * AL 1 
A8--<  A8  +  C22)*AL1 
A10  --C234AL1 
A10-  < A10FC24 )*AL1 
A10  =  <  A10FC25)-»AL1 
A 1 0- <  A 1 0+C26 ) * AL1 
A10  ( A1 0+C27 )  #  AL 1 

B=A10»  T2  +  A8 
B=B*T2+A6 
H-B*T2+A4 
B-BFT2+A2 
TF'H1-B4T2  +  1 .0 

PHI  -TFHI»X-(GNU+.5)*(F- 1/2.0) 

COF‘-'COS  <  PHI ) 

S I F*  '•  S I  H  <  F'  H I  ) 

F1-ANU4C0P 

Y1-ANU4SIP 

IF(KOUNT-l)  10.10.11 

10  FSAVE- F 1 

Y  S A VF -  Y 1 
GNU  FNU  +  1 .0 
F.  0  U  N  T  =  2 

GO  TO  9 

11  F  2 ;  F  1 
Y2:  Y  1 

F  1  --FSAVE 

Y  1  =  Y  S  AYE 

IF(X-GO.O)  110.111,111 

110  IF  (ABS(Fl)  -  AB  S ( F  2  )  >  12,12,13 

111  F  < 1 )  -  F 1 
F  (  2  )  -  F  2 

IF(M)  18,22,112 

112  IF(N-l)  22.22.113 

113  Ml  -NFI+1 
FLI  -  0.0 

BO  114  I ; 3 . N 1 
f  L I  ' FL 1  +  1 . 0 

114  F(  I > -CONST « (FHU  +  FLI >4F ( 1-1  )-F( 1-2) 

GO  TO  22 

12  ALF-  F(2)/F2 
GO  TO  14 
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13  /ILF  - F  (  1  >/Fl 

14  F(  1  >-  FI 
F ( 2) -F2 
11  —  3 

LO  TO  16 
END 

♦DECK  BANDS 

SUBROUTINE  BANSOL ( A i D » B . NR . KH , NUKBl.K i NECT r LLL » E i NLGL » F  ) 
COMPLEX  F 
DIMENSION  F(NN.IiK) 

COMPLEX  A  >  B  >  D  >  C  >  E 

DIMENSION  AtLLL.Hli) .  D  ( 1 ) . B (1 > , E ( NL GL » MM > 

LEVEL  2 » A » D » B  » E  >  F 

NNN  LLL/2 

NG=NN*MM 

NBL:  NN 

MG=N0*2 

MQQ"  MQ 

NBAND-MH 

NBLO^NN 

NB  =  0 

GO  TO  150 

100  l!B  =  NB  +  l 

IF  ( NUMBLK-NB- 1 )  103.104.101 
104  MOO -NNN*MH*2 
NBLC- NNN 
GO  TO  101 
103  NBL - NNN 
C  CHANGES 

101  DO  125  N-l.NBL 
NM-MNTN 

B  (  N  )  =  B  (  N  M  ) 

B(NM)=0.0 
DO  125  M:  1 »  MM 
A(N.M)- A(NM.M) 

125  A  (  N  M  »  M  )  =  0 . 0 
C  NO  CHANGES 

IF  <  NUMPLK  -NB)  150.200.150 
150  NQO-MB+11 

IF  (NBLC-NN)  1 09 i 1 1 0 i 1 09 

110  CALL  FILL  <  A » F  » B .  D .  NB .  NN .  NNN .  NBAND . LLL .  NC!<2 r  MOD) 

GO  TO  126 

109  CALL  FILL  ( A , E . B » D . NB . NN , NNN . 

1  NBANDiLLL.NGCl.NaQ) 

126  IF  (NB)  200.100.200 
200  DO  300  N-l.NBL 

IF  ( CABS  <  A ( N . 1 ) )  .LT.  1.0E-30)  GO  TO  999 
225  B  ( N )  -'  6  ( M >  / A  <  N  >  1 ) 

DO  275  L=2»I'.M 
IF  <  A <  N . L  > )  230.275.230 
230  C-A<N.L)/ACM.l) 

I  =  N  +  L  -  1 
J  =  0 

DO  250  K - L  » M M 
J- J  +  l 

250  A(I. J)  = A(I. J>-C*A(N.K) 

B ( I )  B( I >-A(N.L)*B(N) 

A  (  N  .  L  >  -  C 
275  CONTINUE 
300  CONTINUE 

IF  (NUMBLK-NB)  375.405.375 
375  NOZ  (MB-1 )F11 
NBNN-NB*  NN 
DO  127  N= 1 »  NN 
NC  NBNN-NN+N 
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jj  D(NC)=B(N)  •  j 

HO  127  M- 1  . NBAND 
127  F(N.M)  A(N.M) 

CALL  ECUR<F.NCZ.KG, IERR) 

IF  ( I  ERR  .EQ.  1)  GO  TO  888 
GO  TO  100 
400  CONTINUE 
400  IiO  430  il  l  iNBL 
N=NBL+1-M 
IiO  425  R=2»MM 
L  =  I !  +  K  -  1 

425  B(N)-B(N)-A<N»K)#B<L) 

NH-N+HHN 
KNtHN)-B(M) 

450  A  (  N  H  .  H  B  >  -  B  (  N  ) 

NB-MB-1 
NBL-  K'N 

IF  ( NB )  470(500.475 
473  MQL  (MB-l)tll 

CALL  ECRIKF.NGL.HG.IERR) 

IF  ( I  ERR  .EQ.  1)  GO  TO  777 
DO  134  N=1.NN 
HC-KBCNK-NH+N 
B  (  M  )  -  D  <  N  C  ) 

DO  134  K5: 1  iNBAND 
134  A<N.K)=-F<N.M> 

GO  TO  400 
500  K=0 

DO  600  N6-1.HUMBLK 
IF  <  NB  .EG.  K'UHBLK  >  MBL=NNN 
BO  600  N'-l.NBL 
NM-K+NNN 
K  K  +  l 

600  IKK)-A<NN.NB> 

553  RETURN 
777  K'  =l 

444  FRINT  666. N 

664  FORMAT  (10H  ERROR  NO.. 110) 

STOP 
338  N=2 

GO  TO  444 

99?  PRINT  998. A <  N . 1 ) >  N 

998  FORMAT  (13H  D I  AG  TROUBLE  . 5X . 2G 1 3 . 5 . 15 > 

STOP 
END 

*DECK  FILL 

SUBROUTINE  FILL  ( A . E . B . D . KB . NH . NN . NBD . LLL . NCQ . MQQ ) 

COMPLEX  A.B.D.E 

B I  MENS  I  ON  A < LLL. NBD) . B ( 1 ) . D ( 1 ) .E(NN.NBD) 

LEVEL  2.A.E.B.D 
NBN  NBTNM 

CALL  ECRIKE.NCG.MGG.IER) 

IF  ( IER  .EQ.  1)  GO  TO  999 
DO  1  N -  1 . N  N 
NC= NBH  +  N 
NNN  -NM  +  N 
B  <  NNN ) - D ( NC ) 
no  1  li-l.NBD 

1  A(NNN.M)  -E(N.M) 

RETURN 

999  M=3 
PRINT  666. N 

666  FORMAT  (10H  ERROR  NO.  .13) 

STOP 
END 


»DECK  ECUR 
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SUBROUTINE  ECWR(S,L1 ,L2 , IERR) 

DIMENSION  3(1) 

LEVEL  2»S 
IERR  0 
KLW  I  Nil  LI 

UKITE(Ll)  (S(I).I  liL2> 

1F(UNIT(L1).GE.0.)IERR  =  1 
RETURN 
E  Nil 

*OECN  ECRU 

E'UDROUT  I l<E  ECRD  (  £  i  L 1 ,  L2 . 1 ERD  > 

M MENS  I  ON  S  < 1 ) 

LEVEL  2.S 
IERH  -  0 
REWIND  LI 

READ ( L 1 >  (S(I)iI =1»L2) 

IF  (  (EOF  (LI )  .R'E.O)  •  OR .  (IOCHECCLI  )  .NE,Q)  )  IER0  =  I 

RETURN 

END 
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•  v-.o 


Col. 

1  -  80 


16  -  20 


Format 

20A4 


DATA  SET  NO.  1 
CARD  NO .  1 
Variable  Name 
ID 


1  -  5 

15 

NBAND 

6-10 

15 

IPR 

11-20 

F10.2 

RR 

21  -  30 

F10.5 

ALPHA 

31  -  40 

F10.5 

DM 

41  -  46 

16 

IBLKLEN 

1  -  5 

15 

i 

IPC 

6-10 

15 

IPMAX 

11  -  15 

15 

IDPC 

_ De-script  ion _ 

This  card  contains  alpha-miner i 
title  that  appears  in  tlu 
output  for  job  identifi- 
cat ion 


CARD  NO.  2 


NBAND  equals  either  the  maximum 

difference  between  two  nodes 
or  the  sum  of  NX’S  and  NODR, 
whichever  is  greater. 

IPR=1,  input  data  check  -  echo 
printed 

IPR=2,  echo  not  printed 

Radius  of  the  semicircle  (feet) 

Wave  angle  of  incidence 

Water  depth  of  the  node  on  the 
semicircle 

The  value  of  1BLKLEN  should  not 
exceed  31000 


CARD  NO.  3 


21  -  25 


IBUGA 


1PC=0,  the  elevation  and  velocity 
at  selected  elevation 
stations  are  omitted 

IPC=1,  the  values  are  printed 

IPMAX=0 ,  subroutine  CPUNCH  will 
not  be  called 

IDPC=0,  special  procedure  for 

input  of  average  depths  for 
each  element  is  used 

1DPC=1,  read  depth  at  each  grid 
node  from  cards  or  from 
the  data  files 

IDCH=1,  element  number,  nodal 
point,  and  depth  of  each 
element  and  the  x  and  y 
coordinates  of  each  node 
are  printed 

IBUC.A=1  ,  ELK  (element  stiffness 
matrix  3x3)  are  printed 


.i 

,i 

.N  1 


69 


Col . 

Format 

Variable  Name 

Descript  ion 

26  - 

30 

15 

IBUGC 

IBUGG= 1 ,  SYSK(l.l)  and  SYSK(I,2) 
are  printed 

31  - 

35 

15 

1BUGL 

IBUGL=1,  SYSO(I)  are  printed 

36  - 

40 

15 

ID1SP 

IDISP=1,  dispersion  included 

41  - 

45 

15 

1CIRC 

ICIRC=0,  for  semicircle  grid 
ICIRC=1,  for  circular  grid 

46  - 

50 

15 

NC 

NC=0,  selected  elevation  stations 
(nodes)  will  not  be  punched 

CARD  NO.  4 

1  - 

10 

FI 0 . 2 

TMIN 

The  minimum  wave  period  (seconds] 

1  1  - 

20 

F10. 2 

TDELT 

The  increment  of  wave  period 
(seconds) 

21  - 

25 

15 

IKK 

The  number  of  times  of  the 
increment 

DATA  SET  NO.  2 

CARD  NO.  1 

1  - 

5 

15 

NSP(l) 

Selected  elevation  station 

number  (nodes) ,  starting 
from  NSP ( 1 )  to  NSP(ISPT) 

6  - 

10 

15 

NSP (2) 

Fourteen  values  on  each  card  and 
continuing  to  the  next  card 
as  needed 

66  - 

70 

15 

NSP ( 1 4) 

DATA  SET  NO.  3 

CARD  NO.  1 

1  - 

5 

15 

NSPE(l) 

Selected  elevation  station 

number  (element),  starting 
from  NSPE(l)  to  NSPE(ISPTE) 

6  - 

10 

i5 

NSPE(2) 

Fourteen  values  on  each  card  and 
continuing  to  the  next  card 
as  needed 

70 


Format 


Variable  Name 


Descript  ion 


15 

NSI’E(  1  4) 

DATA  SET  NO.  4 

CARD  NO.  1 

F 1 0 .  V 

X(l) 

x-coordinate  at  each  nodal  point 
starting  from  X ( 1 )  to 
X(NNOD) 

F10.2 

X(2) 

Eight  values  on 
continuing 
as  needed 

each  card  and 
to  the  next  card 

F10.2 

X(S) 

DATA  SET  NO. 

_5 

CARD  NO.  1 

F10.2 

Y(l) 

y-coordinate  at  each  nodal 
starting  from  Y ( 1 )  to 
Y(NNOD) 

point 

F10. 2 

Y(2) 

Eight  values  on 
cent inuing 
as  needed 

each  card  and 
to  the  next  card 

Col 


Format 


Variable  Name 


Desc  r ipt i on 


Blanks 


Blanks 


DATA  SET  SO.  6 
CARD  NO.  1 


NCON (1,1) 


NCOS’ ( 1,2) 


NCON (1,3) 


NCON (2, 1) 


NCON (2, 2) 


Node  number  of  the  first  nodal 
point  of  element  1 

Node  number  of  the  second  nodal 
point  of  element  1 

Node  number  of  the  third  nodal 
point  of  element  1 


vode  number  of  the  first  nodal 
point  of  element  2,  etc. 

Nine  values  on  each  card,  starting 
from  NCON ( 1 , 1 ) ,  NCON( 1 ,2) , 
and  NCON (1,3)  to  NCON (NELE , 1 ) , 
NCON (NELE , 2) ,  and  NCON (NELE, 3) 
respectively,  continuing  to 
the  next  card  as  needed 


NCON (3 , 3) 


DATA  SET  NO.  7* 


F10.2 


F10.2 


CARD  NO.  1 


DD(  1 ) 


DD(2) 


Depth  at  each  nodal  point, 
starting  from  DD(1)  to 
DD(NNOD) 

Eight  values  on  each  card  and 
continuing  to  the  next 
card  as  needed 


F10.2 


DD(  8) 


ote:  If  IDPC«0,  Data  Set  No.  7  is  omitted.  Depth  of  each  node  should  be 

specified  by  user  in  the  program  FINITE. 


49.  To  set  the  correct  values  for  the  parameter  statements  in  program 

FINITE,  the  following  procedure  should  be  followed. 

£.  Sign  on  the  Cybernet  Computer  System 
Telephone  No.  2047  for  300  Baud 
Telephone  No.  2030  for  1200  Baud 

b_.  Type  in 

GET,  FINITE/L'N=CER0M0 
SAVE, FINITE 

c_.  Use  XEDIT  to  set  correct  values  for  the  parameters  in  the 
parameter  statements 
Type  in 
OLD, FINITE 
XEDIT 

d_.  Locate  the  first  parameter  statement 
Type  in 
L/PARAMETER/ 

The  system  will  respond  with 

PARAMETER (NNOD= 3 5 ,NELE=48 ,N0DR=9 ,NCS=7 , ISPT=35 , 1SPTE=48) 

e.  Change  old  values  in  the  parameter  statement  to  desired  values 
Type  in 

C/NNOD=35 ,NELE=48 ,N0DR=9/NN0D=new  value ,NELE=new  value ,N0DR=new 
value/ 

C/NCS=7 , ISPT=35 , ISPTE=48/NCS=new  value , ISPT=new  value  ,  ISPTE=new 
value/ 

f.  Locate  the  second  parameter  statement 
Type  in 

N1 

The  system  will  respond  with 

PARAMETER (NSYSK=1 394, NSYTP=772 ,NSTMP=84) 

g.  Change  the  values  in  parameter  statement  to  desired  values 
Type  in 

C/NSYSK=1 394 ,NSYTP=772/NSYSK=new  val ue ,NSYTP=new  value/ 
C/NSTMP=84/NSTMP=new  value/ 

h.  Replace  program  file  FINITE  by  editted  version 
Type  in 

E,  ,RL 

50.  Establish  an  "AFFILE"  on  the  CYBER  205  System.  An  "AFFILE"  enables 

the  Corps  of  Engineers  (K0E)  user  to  transmit  files  from  KOE  to  the  CYBER  205, 

and  vise  versa. 

£.  Type  in 

GET,AFJCL/UN=CER0MQ 
SAVE, AFJCL 


b.  Use  XEDIT  to  make  necessary  changes 
Type  in 

OLD, AFJCL 
XEDIT 

c.  Replace  205USER  with  actual  user's  ID  of  the  CYBER  205 
Type  in 

C/205USER/USER' S  205  ID/ 

d.  Replace  all  PASSWORDS  with  the  actual  user's  PASSWORD 
Type  in 

C/PASSWORD/ACTUAL  USER'S  PASSWORD/* 

e.  Replace  all  KOEUSER  with  actual  user's  ID  of  the  KOE 
Type  in 

T 

This  command  will  bring  the  pointer  to  the  top  of  the  file 
Type  in 

C/KOEUSER/ACTUAL  KOE  USER'S  ID/* 

_f.  Replace  file  AFJCL  with  the  edited  version 
Type  in 
E, , RL 

_g.  To  execute  file  AFJCL 
Type  in 

GET,ADYJOB/UN=CATHDDC 

BEGIN, ,ADY JOB, AFJCL 

The  system  will  respond  with 

ALL  DONE.  USER  JOB  NAME  IS  COEXXXX 

(Note:  The  status  of  the  job  can  then  be  tracked  by  using  "LINK, ENQUIRE."  When 

the  job  is  complete,  output  will  return  to  the  users  QGET  queue  and  can  be 

accessed  by  using  the  "QLIST"  command.) 

51.  Establish  file  "QLIST",  which  enables  the  KOE  users  to  retrieve 

CYBER  205  dayfiles  and  outputs. 

a.  Type  in 

GET,QLIST/UN=CER0M0 
SAVE, QLIST 

_b.  Use  XEDIT  to  make  necessary  changes 
Type  in 
OLD, QLIST 
XEDIT 

c.  Replace  old  UID  with  new  UID 
Type  in 

C/UID«HASH/UID=XXXX/ 

(XXXX  is  the  user's  index  hash.  It  can  be  obtained  from  "ENQl'IRE(B) . 
XXXX  are  the  first  four  characters  of  the  job  name  that  the  system 
assigned  to  the  user.) 

d_.  Replace  KOEUSER  with  user's  ID 
Type  in 

C/KOEUSER/USER'S  ID/ 


e.  Replace  old  header  with  user's  header 
Type  in 

C/HEADER  1 /USER'S  HEADER/ 

C/HEADER2/USER ' S  HEADER/ 

C/HEADER3/USER '  S  HEADER/ 

f.  Replace  file  QLIST  with  the  edited  version 
Type  in 

E,  ,RL 

52.  Create  an  Update  Library 

a.  Type  in 

GET, F1NUP1 /UN=CER0M0 

BEGIN, RUNE, FINUP1 

The  system  will  respond  with 

15.56.20  SUBMIT  COMPLETE.  JOB  NAME  IS  AJZZXXX 

b.  Check  dayfile 
Type  in 

OLD, DAY  1 
LIST 

53.  Establish  an  Update  Correction  File 

a.  Type  in 

GET , FIND1R/UN=CER0M0 
SAVE, FINDIR 

b.  Use  XEDIT  to  make  necessary  changes  and  replace  the  file  FINDIR 
with  the  edited  version 

54.  Run  program  FINITE  with  Update  Correction  File  FINDIR  (Figure  29). 

a.  Type  in 

GET , F1NI /UN=CER0M0 
SAVE, FIN I 

b.  Use  XEDIT  to  make  necessary  changes 
Type  in 

OLD , FIN  I 
XEDIT 

c.  Replace  205USER  with  the  CYBER  205  user's  ID 
Type  in 

C/205USER/CYBER  205  USER'S  ID/ 

d.  Replace  PASSWORD  with  user's  PASSWORD 
Type  in 

C/PASSWORD/USER’S  PASSWORD/ 

e.  Replace  all  KOEl'SER  with  user's  KOE  USER  ID 
Type  in 

C/ KOEl'SER /USER '  S  KOE  USER  ID/* 

f.  Change  data  file  name  (optional)* 

Type  in 

C/FINDAT1 /NEW  FILE  NAME/ 

C/FINDAT/NKW  FILE  NAME/ 

(*See  Note,  page  83) 
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OLD , FINDIR 
/LIST 

*IDENT  PARAIS 
*D  HA  IN. 2 

PROGRAM  F I NI T  (  T  AF'EI  TAPE  4 ,  T  ARE  8  =TAPF8  ,  OUTRUT  ,  TAPE*  =  OUTFUT  » 

*D  MAIN. 3, HAIM. 4 

1TAPE9-  T  AFE  9  >  T  ARE  1 1  T ARE  1 1 f T ARE  1 2 =TAPE12, 

2 TARE  13  TAFE13.TAPE14 ■ TAPE14 ,TAFE15 =TARE15) 

*D  HAIM. 6 

CCC  SAMPLE  PROBLEM 
»D  MAIM. 8 

C  tt  ************  t*****************t*******»**************** 

*D  HAIM. 14 

» i  t* »***«* »***«***» ^ »«:*»*** ****************************** ********* 

*D  MAIM. 13 

C  NEC1T  =■  TOTAL  NO.  OF  SIH.  LIN.  EOS.  ■-  NNOD  +  NCS 

C  IKK  TOTAL  MO.  OF  CASAS  TO  BE  COMPUTED 

C  J  SR  T -  TOTAL  NO.  OF  SFLECTD  ELEVATION  STATION  WITH  RESPECT  TO  NODES 

C  ISRTE  -  TOTAL  NO.  OF  SELECTED  ELEVATION  STATION  WITH  RESPECT  TO  ELEMENTS 

C  XK  -  WAVE  PERIOD 

C  SYSQ.STHP.AT.AF  HAVE  DIMENSION  OF  NEOT. 

C  X/Y.NOD  HAVE  DIMENSION  OF  NNOD . 

C  XJ.XY.XH.HD  HAVE  DIMENSION  OF  NCS . 

C  ANGL.SN  HAVE  DIMENSION  OF  MODR. 

C  DDrV.U.YC.NCON  HAVE  DIMENSION  OF  NEIE. 

C  NR  AND  2  HAVE  THE  DIEMHSION  OF  THE  NO.  OF  THE  NODAL  POINTS  IN  EACH  CELL 

C  MSP  HAS  THE  DIMENSION  OF  IPST 

C  NSPE  HAS  THE  DIMENSION  OF  IPSTE 

C  IPER  HAS  THE  DIMENSION  OF  5 

C  XK  HAS  THE  DIMENSION  OF  IKK 

C  AP  HAS  THE  DIMENSION  OF  ISRT+1 

C  SYSK  HAS  THE  DIMENSION  OF  NLGLXNB ANDX2+50 

C  SYTP  AND  FTEMP  HAVE  DIMENSION  OF  NL GL XNBAND+ 100 

C  MUNBLK  -NE0T*N6 AND/ I BLKLEN  +  1 

C  NLG  -  NEOT/NUMBLK 

C  MLGL  -  MEOT- ( MUNBLK  -  1)*NLG 

C  NSTMP  NLG+NLGL 

*D  MAIN. 18. MAIN. 21 
C 

PARAMETER (NNOD- 30. NELE  48, NODR  9. NCS *7. ISPT  =35,  ISF  TE  =  4B) 

PARAMETER! NS YSK  1 3 9 4 , HSYTP =772 , NSTMp=84 ) 

DIMENSION  X ( NNOD ) . Y(NNOP) -MOD (NNOD) , XJ(NCS) , 

»XY(NCS) . XH(MCS) . DH ( NCS ) . AHGL ( NODR ) . SM ( NODR > . NSP < 100 ) « I  PER ( 5 ) 

DIMENSION  I D  <  2  0  > »NSPE< 100) 

DIMENSION  Z ( 3 ) 

»D  MAIN. 22, MAIN. 24 

COMMON /L 1 /SYSK  <  NS YSK ) 

COMMON /L2/SY SO (NSTMP) 

C0M;10M/L3/STMR  <  NSTMP ) 

COMHQN/L  4/DD ( NELE ) 

C0MM0N/L5/AT  (  NNOD  +  NCS) 

C0MII0M/L6/AF  (  NNOD  +  NCS  ) 

C0MH0N/L7/SYTF (NSY1P) 

C0MM0U/L9/U( NELE ) rV(NELE) , YC(NELE) 

COMMON/L9/MCOH(MELE,  3  ) 

DIMENSION  FTEMP(K'SYTP) 

»D  HAIM. 25, MAIN. 23 

DIMENSION  VJ( 100) ,VY( 100) 

DIMENSION  MR ( 3 ) 

DIMENSION  AF(IOO) 

*D  MAIN. 32 
»D  MAIN. 32, MAIN. 37 

Figure  29.  Program  FINDIR,  update  correction  file  for  program  FINITE 
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READ  (8,2)  ID 
*0  HAIM. 41 

RE AH <8. 54 35 )  KB AND , I PR , Rft , ALPHA . DM > I BLKLEN 
5435  FORMAT (2I5,F10.2«F10.5,F10.2» 116) 

*D  HAIM. 4 2 
*1  MAIN. 46 

URI  TE  (  6  r  403  )  NUliBLK 

403  FORMAT  < IX, ' NUMBLX  '  ,  15) 

*D  MAIN. 48 
01=1. 

RE AIi  (  8 , 305)  IRC  ,  I  PM  AX  ,  IDPC,  IDCH,  I  BUGA ,  I  BUGG  ,  IBUGl.i  JDISP, 

1ICIRC.MC 

305  FORMAT (1015) 

IF( ICIRC.ME. 1 )G0  TO  3010 

C  FOR  CIRCULAR  REGION  NCS2- < NCS  + 1 )  72  -  TOTAL  NO.  OF  COEFF .  IN  RAD.  DOMAIN 
NCS2  :  ( MCS  +  1 )/2 
NCK2-  2  4-NCS2 

IF  <  NCK2 . NE . NCS )  GO  TO  3010 
UR  I TE  <  6  >  377 ) 

877  FORMAT  < IX, 33HNCS  MUST  BE  ODD  FOR  CIRCULAR  GRID) 

STOP 

3010  IFUPR.EQ.OIGO  TO  876 

876  COMTIMUE 
*D  MAIN. 54 

IF(NLGGtNBAND.LE.NSYSK)  GO  TO  5432 
*D  MAIM. 58 

NNELE  NELE 
NNNOD" NNOD 
NNODR- NODR 
NNCS' NCS 

URITE(6>12)  NNELE, NNNOD, NNODR, NNCS rNEOT, ALPHA, DA2 , DA4 r Rk r D1 , DM 
*1  MAIN. 72 

URITE(6,2010) 

*D  HAIM. 73 
*D  MAIN. 75 

URITE<6,2014)  NUMBLK , NLG , NLGL , I BLKLEN , IPR , NBAND 
IF(MLG.GE.MBAND)  GO  TO  878 
UR  I TE ( 6  >  87?  > 

879  FORMAT ( 1 X , 30HNLG  MUST  BE  GREATER  THAN  NBAND) 

STOP 

878  CONTINUE 
*D  HA  IN. 78 

READ (8, 300)  THIN , TDELT  ,  IKK 
*0  HA  IN. 30 

UR  I TE( 6, 300) THIN, TDELT, IKK 
*D  HAIM. 81 

300  FORMAT (2F10.2, 15) 

*1  MAIN. 92 

IF(IKK.EQ.l)  GO  TO  321 
*1  HAIM. 86 

321  CONTINUE 
*D  HAIM. 37, HAIM. 89 

IF(ICIRC.EQ.l)  DA2  -2.*PAI/N0DR 
IF(ICIRC.EG.l)  DA4-DA2/2. 

*D  MAIM. 90, MAIM. 93 

IF(IPC.EO.O)  GO  TO  3000 
RE  AIK  8 , 306)  (  NSP(  J)  ,  J  -1  ,  ISPT  ) 

READ (8, 306)  (NSPE(J),J  l.ISPTE) 

C  LEVEL  2,U,V,YC 

C  LEVEL  2 , NCON , DD , AF , AT 

»D  HAIM. 94 

306  FORMAT ( 1415) 

*1  HAIM. 100 

T  =  XK ( IK ) 

G  32.2 

IF(IDISP.EO.O)  GO  TO  77 
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78  TM=C*WK*TANH(WK*DM) 

Tl  SORT  (  TM) 

T2  =  2.*F'AI/T1 
T3-XK! IK1-T2 
T4  -  AB  S  ‘  T  3  ) 

IF  1T4.LE . 0 .01  >  GO  TO  74 

DK  -2.  IFAITT3/  O.'K  !  Ik  >*XK<  I K  >  tSQR'T !  G*  I'M)  ) 

UK  -UK  +  DK 
GO  TO  78 

76  ).H  .  54!  1  .+2.*UKtDM/SINH<  2.  *WK*DI'.>  ) 

01  ( >:  M  *  4.1'fAItFftI  )/<T»T*G*UK*UK) 

UR  I TE 1 e  >  74 )  01 

74  FORMAT!//.  2SX , 44HDI SF ERS I ON  INCLUDED.  01  REPLACED  BY  C*CG/G  ‘.F10. 
12) 

77  CONTINUE 
*1  MAIN. 101 

73  FORMAT!//. 25X.24HUAVE  PERIOD  IN  SECONDS  - . F 10 . 2 r 24X . 9HUAVE  NO  =.E 
112.0.//) 

*D  MAIN. 102. MAIN. 103 
*0  HAIM. 106 

l.NELE.NNOO.NODR.NCS . NEQT . NBAND . NCS2 . DD . I DPC . DM . I OCH . IDISF) 

*D  MAIM.  112 

1 . HLG . I  OF  C . I  BUG A. DM . IDISP. T) 

*0  MAIM. 118 

1 ,NELE .NNOO, NODR. NCS .NEOT. NBAND. NCS2.ND0R1 . VJ.VY.NL GL . I BUGG  > I  CIRC) 
*D  MAIN. 124 

1 .NELE. HNOD, NODR. NCS. NECT. NBAND .NCS2. NODR  1. IBUGL. ICIRC) 

*0  MAIN. 131 .MAIN. 138 
»D  MAIN. 145 

IF!  ICIRC. NE. 1)  AF (I)-AE(I)/2. 

*  I  MAIM. 191 

C  URITE(6.37)  <  I , AF ! I  )  .  At! I ) . I  l.NNOD) 

C  37  FORMAT ! 4 ! IS >  2F7 . 2 ) ) 

*D  MAIN. 197 

240  2 ( L  )  -  AF ! N  ) 

*1  MAIN. 21 2 
JF  -1 
*1  MAIN.  213 

JF'  JFP2 

*0  MAIN. 218. MAIM. 219 

WRITE (4 , 2S )  Jr AF ! J) .AT ! J) . J1 .AF! J1 > . AT ( J1 ) 

501  JFP-JP+1 

A  F' !  J  F' )  AF(J) 

IF ( II  ,GT  .  ISPT  >  GO  TO  S03 
AF' !  JF  F' )  -AF  !  J1  ) 

*0  MAIM. 221 

503  WRITE!*. 24)  J . AF ! J ) . AT ! J ) 

*1  MAIM. 222 

IF  <  M C . E Q  >  0  >  DO  TO  103 
UR  I  TE ! 6 . 101 ) 

101  FORMAT ( 4 H  CCC) 

103  URITET4.201  > 

201  FORMAT'//) 

WRITE  <  4 . 200 )  ( AF ( I  I  )  .  I  I  -  1 . I SF  T  > 

200  FORMAT! 13F4.2) 

*D  MAIN. 323 

502  CONTINUE 

»0  HAIM. 22 4, NAIM. 234 
*D  MAIN . 240. MAIN . 241 

URITE!A.2S)  I . AF ' I ) . AT i I ) . I  I . AF  < I  I ) , AT ! I  I ) 

#0  MAIM. 213 

4 1 S  WRITE! 4.21)  I . AF ! 1 > . A T ! I ) 

*0  MAIM. 24S. MAI. I.2S5 

222  FORMAT ! 3 ( IS. 1 F.3E12. 5. IX ) ) 

*0  MAIN. 257 

IF(NC.EO.O)  GO  TO  104 
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WRITE  (6,102) 

102  FORMAT ( 1H  CCL) 

101  IF(IPC.EQ.O)  GO  TO  99 
*D  MAIN. 278 .MAIN. 280 
*£i  MAIN. 282.  MAIN. 285 

24  FORMAT < IX. 15, 2E 12.5) 

25  FORMAT  <  1X,2(  15.  £12.5.2):.  El  2.5)  ) 

C  IF  <  IF'MAX  ■ NE . 0 )  CALL  CRUNCH  <  NODF  T  ,  IF  MAX .  AF  ,  U.  V  r  NNOD .  NEl.E ,  MFCT ,  XK ,  IK , 

C  1NSP. ISRT. I  PER, ID. IPC.MSRE. ISFTE.IKT) 

*D  INPT.3 

1 . NELE .  NNQD »N0DR,NCS.NEQTiN6AND,NCS2.DP,  I  DF  C .  DM .  I DCH ,  IDI SP ) 

*1  I  MPT.  8 

C  LEVEL  2.NC0N.DD 

*D  INPT.IA.IMPT.37 

READ  <  4 , 200  )  < X < I ) , I -1 , KNOD ) 

RE AD ( 4 , 200  >  ( Y ( I ) , 1=1 , MMOD ) 

200  FORMAT  <  8F 1 0 . 2 ) 

C  DO  1000  1=1 .MMOD 

C  X(I>=X(I)*.lltO 

Cl  000  Y<  I) -YU  >*.11 10 
C  DO  647  1*1055. NNOD 

C  II--I-N0DBI  +  1 

C  IF(I.LT.llOO)  AAD- ( I -  1 055 ) *DA2  + . G*DA2 

C  IF  <  I  .  GE  .  1 100  )  AAD  -*  ( 1 1  —  1 )  *DA2 

C  RS 1 1 60000  . 

C  IF(I.LT.llOO)  RS= 140000 . 

C  :<<I)  =-RS*COS< AAD+ANGG) 

C  Yd)=-RS*  SIN  (AAD  +  ANGG) 

C  647  CONTINUE 

101  FORMAT  <3(5 X. 315) ) 

READ (4 , 101 )  ( <  NCON  < I » J ) , J  =1 r  3 ) > I  1 .NELE) 

C  DO  83  I- 1, NELE 

C  83  DD<  I  >•- 6.*DD(I  > 

CALL  GDEPTH ( I DPC  >  DD »  NNOD  >  NQDB2 , NELE , I DCH ) 

IF< IDCH.EQ.l )  GO  TO  430 
GO  TO  450 

430  PRINT  316, (I ,X( I) , Y< I ) ,1-  1  .NNOD) 

DO  440  1=1, NELE 

PRINT  420,  I.  <NCON<  I  ,  J)  ,  J=-l  ,3) 

440  CONTINUE 

316  FORMAT (2(110. 2F 10. 2)) 

420  FORMAT (IX, 6( 15. IX >315) ) 

450  IF(IDISP.OT.O)  GO  TO  70 
IF(IDPC.NE.O)  GO  TG  400 
*1  INPT.38 

DD  < I ) - DD ( I ) /DM 
*D  INPT.44.INr-T.48 
70  CONTINUE 
*D  INF' 1.4  9 

416  F0RMAT(4(4X,4I5) ) 

*D  IMPT.50 

IF( IDCH.EQ.  1)  PRINT  416,  < L , ( NCON ( L , J ) , J  1 , 3 )  ,  L = 1 , NELE ) 

*1  DEPTH. 7 
C  LEVEL  2. DP 

*1  DEPTH. 9 

DO  10  I  - 1 >  NELE 
10  DD < I ) -0 . 75 
*D  DEPTH. 10, DEPTH. 40 
*D  DEPTH. 45, DEPTH. 46 
1000  READ (4 ,980)  ( DD ( I ) . I = 1 , NNOD ) 

990  FORMAT (8F10.2) 

IF (IDCH.EQ.l) PRINT  2010,  ( DD ( I ) , I = 1 . NNOD ) 

*D  DEPTH.  17 
*1  BAND. 4 

C  LEVEL  2, NCON 

*D  ASBK.5 
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1 ,NLGG, IDPC, IBUG,DM, IDISP.T) 

*0  ASBK.21 

C  LEVEL  2.SYSK.SYSQ.ST.NC0N.DD 

*D  ASBK.73 

C  IFCL.GE. 14  33. AND. L.LE. 1571 >D9- 0.5 

*B  ASBK.80 

CALL  ELMK<  XI . X2. X3 . Y1 »  Y2  >  Y3.UK .D9>  DIN.EL  N. AREA.l.  .KELE3.DH. IDISP. T ) 
*D  ASBK.84 
1030  DO  20  1-1.3 
*1  ASBK.104 

LA  IQ  +  9 
*B  ASBK.113 

999  URITE'6.777> 

*D  GSPK.  3 

1 . MELE . NNOD  >  MODE >  KCS  > MEQT  » NB  AND » NCS2 r  NODR  1 . V J . V Y  >  MLGG > I  BUG. 1C  IRC) 

*B  GSPK. 14 
*D  GSPK. 12 

1 >  ANGL  (  NODR  ) 

*1  GSPK.22 

CALL  BESSY(UR.O. .MCS.VJ.VY) 

NDUM=NCS 

IF< ICIRC.EG. 1)  NDUM-NCS2 
IF ( I  BUG  .  EQ . 0 )  GO  TO  13 
WRITE<  6.21 ) 

21  FORMAT (//. IX. 'XJ- ' . 10X. 'XY  •  '  ) 

*D  GSPK.23 

13  DO  12  I-l.NDUM 
*D  GSPK.23iGSFK.26 
*1  GSPK.31 

IF ( I  BUG . EQ • 0 )  GO  TO  14 
UR  I TE ( 6 . 22 ) 

22  FORMAT <//. IX. ' XH- ' , 10X. '  DH=  ' ) 

*B  GSPK.32 

14  DO  13  I-l.MDUM 
*D  GSPK. 43 

20  FORMAT  ( 9  H  GSPK  BY--.4E15.8) 

*D  GSPK. 30 

DO  38  1=1. HCS 
II  I-  1/2  +  1 
*1  GSPK. 55 
ICE  =1 

IF  < ICIRC . EQ . 1 >  ICE  =  I 1 1 
IF< ICIRC.EG. 1 >  C  X =  2 . *  C  X 
*D  GSPK . 56 

SYSIsdl  ,  1  )=-SYSI<<  II.  1  >+CXtRKHF*XH(  ICE)*DH(  ICE) 

*1  GSPK. 60 
CT=  1  . 

IF< ICIRC.EQ. 1 )  C  T  =  2 . 

ND -M0DR1 

I F  <  ICIRC.EQ. 1)  ND  NODR 
»D  GSPK. 61 

DO  40  II. ND 
*1  GSPK. 69 

IF(Il.EQ.O.AND.ICIRC.EG.l)  II  NODR 
*1  GSPK. 71 

IF ( ICIRC .EQ. 1 )  GO  TO  933 
GO  TO  936 
933  J2= J/2 
J1  J2  +  1 
Al-  A N G L < I >*J2 
A2  ANGL  < II )»J2 
936  CONTINUE 
*D  GSPK. 74 

TH'CT  <  HKL  FDH ( 1 ) 

IF< ICIRC.EQ, 1  GO  TO  46 
»I  GSPK. 77 
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IF  (ICIRC  .EQ.  1)  GO  TO  937 

00  TO  938 

937  IF  <  MOD  <  J  >  2  > • EQ • 1 )  GO  TO  44 

TM  HKL*DH< J1 )*<C0S(A2)+C0S< Al> ) 

00  TO  46 

938  CONTINUE 
*1  GSFK.88 

GO  TO  48 

44  TIV-HKLTDHl  J1  > *  1 S IN  W.2 > +S IN <  A1  >  > 

*D  GSPK.91 
*D  GSPK.97 

C  UKIT£<6»  1018)  (SVSK  ( 1 1 1 )  r  SVSK<  I  >  2)  r  I-=l  i  NLGG ) 

*D  GSPK.98 

C1018  FORMAT  < 1 X » 4G 1 4 . 5 ) 

*D  GSFK.100 
2000  WRITE!  8*999) 

*0  GSPK  . 103 

1551  WRI TE<  6 1 1520 )  NLLU 
»D  LOAD. 3 

1 t MELS.NNOD. NODR t NC3 , NEOT . NBAND » NCS2 f H0DR1 » IBUG » ICIRC) 

»D  LOAD. 9 

DIMENSION  SH(NODR) i SYSO(NEQT) >DH(NCS) >  X J  <  HCS ) »AHGL(NODR) 
*D  LOAD. 10 
*1  LOAD. 15 

NDUM- N0DR1 

IF  < ICIRC . EQ . 1 )  NDUM  NODR 
*D  LOAD. 18 

DO  12  I -  1 1 NDUM 
*1  LOAD. 23 

IF ( ICIRC . EC. 1 )  DUH 2  (O.OiO.O) 

*1  LOAD. 31 

IF(ICIRC.EQ.l)  GO  TO  52 
*1  LOAD. 33 

52  CONTINUE 
*1  LOAD. 34 

IF  ( II .  EQ.O.  AND.  ICIRC.  EQ.  1 )  II --NODR 
*1  LOAD. 19 

IF  < ICIRC . EQ . 1 )  II  1/2 
1 J-II+1 
»I  LOAD. 51 

IF ( ICIRC .EQ . 1 >  GO  TO  555 
GO  TO  558 

555  IF(I.EQ.l)  GO  TO  557 

IF (MODI  1 .2) . EC . 1 )  00  TO  18 
Th  CM»XJ< IJ)#DHCI J)*COS< AI) 

GO  TO  20 

557  TH» XJ< 1 )»DH<1 ) 

GO  TO  20 

18  TK= CM*XJ(IJ)*DH< IJ)»SIN(AI> 

GO  TO  20 

556  CONTINUE 
»I  LOAD. 32 

20  CONTINUE 
*D  ELMT.3 

INELE3.DM. IDISF.T) 

*1  ELMT.25 

C  URITE(6.500)  X 1 1  X  2  1  X  3  >  Y 1 » Y  2  >  Y  3 

*1  ELNT.30 

XK3 • WK/SQRT (D) 

IF<  IDISP.HE.  I )  GO  1(1  300 
XK1  -XK1*S0RT(DM) 

PI  -3.1416 

150  TM-G*XKl*TANH(XKl*in 
Tl-  SQRT(Trl) 

T2 ‘2.4F  I/T1 
T3=T-T2 
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T4*ABS(T-T2) 

If <T4.LE.  .01)  GO  TO  200 
DK"-2.*PI*T3/<T*T*SQRT<G*D)  > 

XK1  -XK1  +  DK 
GO  TO  100 
200  CONTINUE 

XN-  .  3*  <  1 .  +2.  *XKi  *D/SINH<2.*XK1*D) ) 

XD*  <XN*4  .*PUPI)/<T»T*G*XK1*XK1) 

GO  TO  -100 
300  XD-D 
400  CONTINUE 
tD  ELHT .31 iELMT . 33 

UK A-DIN*XK1«XK1* AREA/12. 

SAU  -D 
D=XD 
*D  ELHT. 43 
D  =  SAV 
*1  BEST. 11 
FN=  N 
*D  BESY.80 

G  A  M  3  = -GAM3TGAM2 
*D  BESJ.38 

5  f  (I ) 'CONST* <FLUFNU)*F<  I+l)-F<I+2> 

*D  BESJ.76 

21  F<L)*C0NST*FRAC*F(L-l)-F<L-2> 

»D  BES J .202 

114  F<I  )=-CONST*(FNU  +  FLI  >*F<  1-1  >-F<I-2> 

*D  GAMA.03fGAHA.36 

XFT  =EXP(TEH)/>:FT 
AMS’ XFT 
GO  TO  6000 
4000  SUM- 0.0 
*1  BANDS. 5 

DIMENSION  QUOOfIOO)  fRUOOf  100)  fMLUOO)  fHLOOO) 
*0  BANOS. 7 
*D  BANDS. 44 

IF(REAL<A<N.L) ) )  230f273f230 
*D  BANDS. 62 

F  <  N  f  M  )  *  A  <  N  f  M  ) 

127  CONTINUE 
*D  BANDS. 74 

A<NMfNB)-'B<N) 

450  CONTINUE 
*D  BANDS. 85 

A  <  N  f  H  )  =  F  <  N  f  M  ) 

136  CONTINUE 
*0  BANDS. 93 

D<K)=A<NMfNB> 

600  CONTINUE 
tD  BANDS. 96 

444  URITE<6f666)N 
*D  BANDS. 101 » BANDS . 102 

999  URITE<6f998)  A<Nf1)fNfNB 

998  FORMAT  < 13H  DI AG  TROUBLE  f 3X . 2G 1 3 . 5 . 2 15 ) 

*0  FILL. 17 

URITE<6f666)N 
*D  FILL. 5 

DIMENSION  Q< IOOf 100) fR< 100- 100) fNL(IOO) fML< 100) 
*D  ECUR. 4 
»D  ECUR . 7 i ECUR . 9 

BUFFER  0UT<L1 fO) <S<1 > iS<L2)  ) 

IF  <  UNI f <L1 ) ) 10  f20 i 30 
30  WRITE(6f40> 

40  FORMAT! 1 X , 'PARITY  ERROR') 

20  IERR'1 
10  RETURN 
*D  ECRD.A 
*D  ECRD.7.ECRD.9 

BUFFER  IN<L1 fO) <S< 1 > fS<L2) ) 

IF  <  UNIT  <L1))10f20f30 
30  WRITE(6f40) 

40  FORMATIIXf 'PARITY  ERROR') 

20  IERD«1 
10  RETURN 
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(*Note:  File  FINDAT1  contains  Data  Set  No.  1  through  Data  Set  No.  3.  File 
FINDAT  contains  Data  Set  No.  4  through  Data  Set  No.  7.  File  FINDAT  can  be  ob¬ 
tained  at  the  completion  of  the  run  of  the  plotting  program  FNGRID.  If  that 
is  the  case,  the  new  file  name  replacing  FINDAT  should  correspond  with  the 
file  name  that  was  used  to  save  TAPE10  in  the  procedure  file  PLOTF.) 

g.  Replace  file  FINI  with  the  edited  version 
Type  in 

E ,  ,RL 

h.  To  execute  file  FINI 
Type  in 

SUBMIT, FINI 

i.  To  check  job  status 
Type  in 

LINK .ENQUIRE , PD=MODD , PT=HHMM 

(PD=MODD,  where  MO  is  the  month  from  01  to  12,  and  DD  is  the 
date  from  01  to  31.  PT=KHMM,  where  HH  is  the  hour  from 
01  to  24,  and  MM  is  the  minute  from  01  to  60.  Example: 

To  generate  status  information  for  a  job  received  since 
9:00  a.m.  on  31  May,  type  in 

LINK, ENQUIRE, PD=0531 ,PT=0900 

The  last  three  characters  of  the  run  JOB  NAME  will  not  stay 
the  same  as  in  initiated  job  names.  The  system  assigns 
different  names  at  different  stages  that  transmit  a  job 
from  KOE  to  CYBER  205,  or  vice  versa.  If  the  job  output 
is  available  at  KOE,  then  continue  to  the  next  step.) 

j_.  To  execute  "QLIST" 

Type  in 
OLD, QLIST 

BEGIN, , QLIST, XXX, LSV.DSV 

(where  XXX  are  the  last  three  characters  of  the  job  name  that 

appeared  in  the  statement  "JOB  NAME  IS  IN  USER  QGET  QUEUE." 
LSV  causes  the  output  to  be  saved  with  file  name  HASHXXX. 
DSV  causes  the  dayfile  to  be  saved  with  file  name  DAYXXX. 

k.  To  list  output  and  dayfile  on  the  user's  terminal 
Type  in 

OLD, FILE  NAME 
LIST 

l.  To  direct  output  and  dayfile  to  a  fast  printer  at  COPE  terminal 
Type  in 

GET,LFN=FILE  NAME 

ROUTE, LFN,DC=PR,UN=C0PE  USER’S  ID 

(where  LFN  is  the  local  file  name,  and  FILE  NAME  is  the 
permanent  file  name) 


55.  In  this  sample  application  to  demonstrate  the  execution  of  program 


FINITE  on  the  CYBER  205  COMPUTER  SYSTEM,  the  values  of  the  parameters  in  the 
parameter  statements  stay  unchanged. 

a.  Obtain  the  source  file  FINITE 
Type  in 

GET ,FIN1TE/UN=CER0M0 
/SAVE, FINITE 

b.  Establish  AFFILE  "AF205"  on  CYBER  205 
1.  To  execute  file  AFJCL 


GET  »  AFJCL/UN-CEROMO 
/SAVE* AFJCL 
/OLD. AFJCL 
/XEDIT 

XEPIT  3.1.00 
??  C/20GUSER/052010/ 
USER(U-0S2010»PA-  PASSWORD) 

??  C/ PASSWORD/ IODINE/* 

USER (U-052010» PA- IODINE) 

--EOR — 

USER ( KOEUSER . IODINE . KOE ) 

END  OF  FILE 
??  T 

??  C/K0EUSER/CER0M8/* 

CHARGE (CER0EGC.CER0M8) 

--EOR-- 

USER(CER0N8* IODINE »KOE> 

CHARGE (CER0EGC.CER0M8) 

END  OF  FILE 
??  Ef, RL 

AFJCL  REPLACED 
AFJCL  IS  A  LOCAL  FILE 
AEB  .  0 , 2S7UNTS . 

/GET , ADY JOB/UN=CATHDDC 
/BEGIN. . ADYJOB. AFJCL 
ALL  DONE.  USER  JOB  NAME  IS  C0E429 
$REOERT . 


2.  Use  "LINK .ENQUIRE"  to  check  job  status 


LINK  .EN0U1FE  .  JN-CCt  429 

STATUS  F  OF;  CENOKB  /KOE 


GATE 

TINE 

USEE 

S  Y  S  T  E  H 

FILE 

FILE 

•  ,L 

n  n  o  v 

HHNM 

JOB  NAME 

JOE  il.TNE 

TYF  £ 

STATUS 

0022 

1C  2  0 

COE  42? 

A  J2Z302 

NO 

TOUTING 

INITIATE G  TO  KAA 

0222 

1  244 

COE  429 

SENI'LOT 

TO 

AF.F.IUEG 

AT  AGY 

• 

0222 

12  02 

CCE429 

AUU.V  00 

NO 

AFRI'.'EG 

AT  K  AA 

0222 

1242 

C0E429 

AUUALOO 

TO 

AUUALCO 

LINKING  TO  AI'Y 

0222 

LINK 

1244 

C  G . :  f 

CCE429 
lETE  . 

A JZZLCT 

WT 

DUTF'UT 

AVAILAHE  FT  KOE/T 

84 


3.  Use  "QLIST"  to  access  the  output 


OLD. OUST 
/BEGIN  » .QLIST >  LOT 

11.44.47  RESOURCE  <  JCAT  P4.TL=200> 

11.44.49  CHARGE  »  CEROEGC .  CEROHB 

11.44.49  ALL  DOME 

11.44.49  T V < 10  +  . 

11.41.49  PURGE* AE205) 

11.44.50  ALL  HONE 
11.44.50  TV . 4+ , 

11.44.50  COPT . INPUT . AF 205 • 

11.41.51  9  WORDS  OF  FILE  INPUT  COPIED  TO  FILE  AF205 


11.44.51  ALL  DOME 
11.44.51  DEFINE. AF205. 

11.44.51  EXISTING  LOCAL  FILE  MADE  PERMANENT 
11.44.51  ALL  DONE 

11.44.51  CHARGE. CEROEGC  .  CER0I18  3.982 

11.44.51  SYSTEM  BILLING  UNITS  (SBU)  3.982 

11.44.51  USER  CPU  TIME  (SECS)  .194 

11.44.51  SYSTEM  CPU  TIME  (SECS)  .493 

11.44.51  USER  MEMORY  USAGE  ( PAGEBSECS )  18.731 

11.44.51  USER  AVERAGE  WORKING  SET  SIZE  (PAGES)  94 

11.44.51  NUMBER  OF  VIRTUAL  SYSTEM  REQUESTS  235 

11.44.51  HUMBER  OF  SMALL  PAGE  FAULTS  73 

11.44.51  NUMBER  OF  DISK  I/O  REQUESTS  11 

11.44.51  NUMBER  OF  DISK  SECTORS  TRANSFERRED  11 

11.41.51  MCOMPLETEt* 

END  OF  FILE 

/ 


SBUS 


c.  Establish  the  file  "QLIST" 


GET .OLIST/UN-CEROMO 


/SAVE. QLIST 
/OLD. QLIST 
/XEBIT 

XEDIT  3.1.00 
7?  C/HASH/AJZZ/ 

.FROC.CLIST .FID.LSV  N/Y.DSV 

IpROC*.  CLIST^FID^LSV^  N/Y»DSV*N/Y.OSP'H/Y»UID*AJZZ» UN-CER0M8. AUH 
77  C/HEADER 1/LUCYCHOU/ 

HI =  LUCYCHOU .H2  =  HEADER2 .H3-HEADER3. 


H/Y.DSF-  N/Y.UID -AJZZ.UH  USERNUM.AUH  --'•HENDDCf 

CHENDDC. 


77  C/HEADER2/CERC/ 

HI ^  LUCYCHOU. H2  =  CERC . H3 -HE ADER3 . 


77  C/HEADER3/UE3/ 

Hl-L UCY CHOU. H2CERC.H3;  WES. 


7?  E..RL 

OLIST  REPLACED 
OUST  IS  A  LOCAL  FILE 
AEB  .  0.251UNTS. 

/ 


I 


d.  Create  an  Update  Library 

GET » FINUP1/UN=CER0M0 
/SAVE.FINUP1 
/BEGIN »RUNF  » FINUP1 

15.56.41.  SUBMIT  COMPLETE.  JOBNAHE  IS  AJZZFOT 
♦REVERT . CCL 

/ENQUIRE. JN=F0T 
AJZZFOT  NOT  FOUND. 

/OLD »  DAY 1 
/LIST 

15.56.41 . FI NI JOB. CM320000 . P6 . 

15. 56. 41.  ABF  .  INPUT  .  0 . 002KI0DB . 01 . 

15.56.41 . USER.CER0M8. >KOE . 

15. 56*41.  ABG  .  F'6. 

15.56.42. CHARGE.CER0EGC.CER0M8, 

15. 56. 43.  GET. FINITE . 

15.56.44. UPDATE(I-FINITE»N::LIBLC) 

15.56.44 . UPDATE  CREATION  RUN 

15. 56. 44.  CREATING  NEW  PROGRAM  LIBRARY 

15.56.46.  UPDATE  COMPLETE. 

15.56.46.SAVE.LIBLC -LIBLC. 

15.56.47. DAYFILE.DAY1. 

/ 

e.  Establish  an  Update  Correction  File  "FINDIR".  (In  the  sample 

run,  file  FINDIR  is  used  without  change.) 

Type  in 

GET,FINDIR/UN=CEROMO 
/SAVE .FINDIR 
/ 


f.  To  execute  program  FINITE  with  Update  Correction  file  FINDIR 
1.  Execute  file  FINI 

GET  ,  F  I M I  / U N  - C E F; 0 M 0 
/SAVE. FINI 
/OLD » F I  III 
/  >:  E  D I T 

/EDIT  3.1.00 
??  C/205USER/052010/ 

USER ( U  052010: PA- PASSWORD ) AD Y 
??  C /PASSWORD/ IODINE/ 

USER (U -052 01  Or  PA  IODINE) AD Y 
??  C / 1 . 0 E U S E R / C E R 0 1 1 8 / T 
CHARGE (CEROEGC .CEROMS) 

GETKOE  >  C0HPF/DD--C6  . UK  =  CEROMB  . 

LINK.  REPLACE  (LIST  LIST205/UI!  CEROHS  .  FM  -K0E  :  DD  C6  .  AF  AK205) 
LINK  ,  REPLACE  (  F  H  M  A  P ' :  M  A  P  2  0  5  /  U  M " :  C  E  R  0  ri  3  .  FM=  KOE  .  DD-  Co.  AF  AF205  ) 
LINK. REPLACE (LIST  LIST205/UN  CEROHS . FM  -ROE ; DD  C6 . AF  Ar  205) 
LINK  .  REPLACE  (  F l-MiAF-1  MAP205/UI!  CEROHS  .  FH  -KOI  .  DD-CS »  AF  A!-  205  > 
END  OF  FILE 
'i”1’  E..RL 

FINI  REPLACED 

FINI  IS  A  LOCAL  FILE 

AEB  .  0.255UMTS. 

/SUBMIT .FINI 

14.24.02.  SUBMIT  COMPLETE. 


JOBNAMF  JR  AJZZWQA 


2.  Check  job  status 

LINK.ENQUIR£*PD=0531»PT=0100 

STATUS  FOR  CER0H8  /KQE 


- 

DATE 

TIME 

USER 

SYSTEM 

FILE 

FILE 

r 

MHDD 

HHMH 

JOB  NAME 

JOB  MAHE 

TYPE 

STATUS 

r 

0531 

1421 

C0E553 

AJZZUQA 

NO 

ROUTING  INITIATED  TO  KAA 

L 

0531 

1433 

C0E555 

UUCYG1 J 

TO 

ARRIVED  AT  ADY 

r* 

0531 

1433 

C0E555 

AUUAGIH 

MO 

ARRIVED  AT  KAA 

0531 

1434 

C0E555 

AUUAGIH 

TO 

AUUAGIH  LINKING  TO  ADY 

0531 

1439 

C0E555 

AJZZGI J 

UT 

OUTPUT  AVAILABLE  AT  KOE/T 

0 

LINK 

COMPLETE. 

3.  To  retrieve  job  output  and  save  it  as  a  permanent  file 
Type  in 

OLD, QL I ST 

/BEGIN, ,QLIST,GIJ,LSV,DSV 

4.  To  route  the  output  to  the  remote  COPE  terminal 
Type  in 

GET ,OUTLUCY=AJZZGI J 
/ROUTE , OUTLUCY , DC=PR , UN=CERORC 
The  system  will  respond  with 
ROUTE  COMPLETE 

/ 

5.  A  lis  .ing  of  the  sample  application  output  from  the  COPE 
terminal  includes  the  following: 


13.29.04 

13.29.04 

13.29.04 

13.29.06 

13.29.06 

13.29,04 

13.29.07 

13.29.14 

13.29.14 

13.29.14 

13.29.14 

13.29.16 

13.29. 16 

13.29.17 

13.29.17 

13.29. 18 
13.29.13 
13.29.24 
13.29.24 
13.29.24 
13.29.24 
13.29.21 
13.29.24 
13.29.24 
13.29.24 
13.29.24 
13.29.21 
13.29.24 

13.29.24 
13.29.23 
13.29.23 
13.29.23 
13.29.23 
13.29.23 
13.29.23 
13.29.23 

13.29.25 
13.29.23 
13.29.23 

13.29.23 

13.29.26 
13.29.26 
13.29.26 
13.29.26 
13.29.26 
13.29.26 
13.29.26- 
13.29.26 

13.29.24 
13.29.26 

13.29.26 

13.29.26 

13.29.27 
13.29.37 
13.29.37 
13.29.41 
13.29.41 
13.29.41 
13.29.41 


*  4  * 

1 4  4  0E  TKOE ( TAKE  4 •  PINDAT) 

4-44 

MFLINKL  T.'.FE4.ST  KAA»DD=C6» 

JCS  ’  4  4  4  44  4  44  *4 9: * 4: 444 444*4:  *  i 
’4444444444  444  4444’  ) 

WAITING  FOR  CONNECTED  STATUS. 

RTFS  -  USER  ( CQEAPPL » ) 

RTFS  -  ATTACH(C0E553/NA> 

RTFS  -  FF  RECUEST  COMPLETE. 

ALL  DOME 

FORTRAN (I  C0I1FF.L  LIST) 

FORTRAN  2.1.3  CYCLE  L592C  BUILT  05/01/81  16121 
COMPILING  FIMIT 
MO  ERRORS 
COMPILING  CF'UMCH 
NO  ERRORS 
COMPILING  INPUT2 
NO  ERRORS 
COMPILING  GDEPIH 
NO  ERRORS 
COMPILING  HAND 
NO  ERRORS 
COMPILING  A3EMBK 
MO  ERRORS 
COMPILING  GSRK 
NO  ERRORS 

FILE  BINARY  EXTENDED i  liEU  LENGTH  ■■■  32 

COMPILING  l  OAD 
MO  ERRORS 
COMPILING  EL  MR 
NO  ERRORS 
COMPILING  BESSY 
NO  ERRORS 
COMPILING  GAMMA 
NO  ERRORS 
COMPILING  BESS J 

NO  errors 

COMPILING  BANS, ML 
NO  ERRORS 
COMPILING  FILL 
NO  ERRORS 
COMPILING  ECUR 
MO  ERRORS 
COMPILING  ECRl) 

MO  ERRORS 

1.914  SECONDS  COMPILATION  TIME 

ALL  DONE 
LOAD i L  FNMMP . 

LOAD  R2.1  CYCLE  L592C 
ALL  DOME 

GO ( T  APE  6  OUTPUT) 

41  STOP  44 
ALL  DONE 

LINK . REPLACE ( LIST  LI  ST  203/UN ;CER0N8  t  FM  KOF  »  DD  C6r  AF  AF205) 

LINK  PHASE  II  VERSION  LLINKOO  BUILT  05/01/81  11.47.59 
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13.39.11  RtFLACE<LIST205/llN  CEROMH. 

13.29.41  L38B2.F6.  KOE  JOB  1.LINK50  3882 
13.29.T.1  WAITING  FOR  CUMNECTFD  STftlUS. 

13.30.03  IKANSFER  ABORTED 

13.30.10  ALL  DOME 

13.30.10  LINR.REPLACEIFHMAP  MAP2G5/UN  MKROMB.FH  KOF.DD  C6.A  y  AF  200) 

13.30.11  LINK  PHASE  II  VERSION  LL1NK50  BUILT  05/01/111  II. 47.5V 

13.30.11  REPl  ACE  (  HAF  205/UM  CER0H8. 

13.30.12  L3EU2.PA.  KOE  JOB  LLINKSO  3LU2 

13.30.18  WA11IMG  FOR  COMUECUD  STATUS. 

13.30,31  TRANSFER  ABORIED 

13. 30. 3 A  ALL  HOME 

13.30.36  I, FI  I  UK (HULL  i ST  KAA.DH-C6. 

13.30.36  JCS-'IEILII  !  6  4  *  4  4  4  4  4  4  4;  4-  4  4  P  *  » 

13.30.36  '  iHitfUit  Z 

13.30.10  WAITING  FOR  CONNECTED  STAIUS. 

13.30.19 

13.30.39  RTFS  -  USER  <  COEA.PF  L  i  ) 

13.30.19 

13.30.19  P1FS  -  PURGE  <  COF  558/NA) 

13.30.19 

13.30.99  F'T  FS  -  PF  REQUEST  COMPLETE . 

13.30.50  ALL  LOME 

13.30.50  MFL 1 44  K  <  HULL  »  ST  KAAiDD=C6i 

13.30.50  JCS- "*4M 4 ****««*f ********** 

13.30.50  •**********»**#»«♦) 

13.30.52  WAITING  FOR  CONNECTED  STATUS. 

13.30.53 

13.30.58  RTFS  -  USER ( COEAPPl. »  ) 

13.30.53 

13.30.58  RTFS  -  PURGE ( CCE557/NA ) 

13.30.53 

13.30.58  RTFS  -  FF  RECUEST  COMPLETE. 

13.31.00  ALL  DOME 

13.31.00  NFL  INK  ( HULL  i  ST  KAA.i  DD-C6. 

13.31.00  JCS-*  *4**  4'I  6*4  **#!***«•****•» 

13.31.00  *  4  **  4  *>  M  >;:*** *****  •  ) 

13.31.01  WAITING  FOR  CONNECTED  STATUS. 

13.31.08 

13.31.08  RTFS  -  USER ( COE AF PL f ) 

13.31.08 

13.31.08  RTFS  -  PURGE ( COF  556/M A ) 

13.31 .08 

13.31.03  RTFS  -  PF  RECUEST  COMPLETE, 

13.31.09  ALL  DOME 


13.31.09  CHARGE .CEROEGC  .CER0M8  71.367  SKUS 

13.31.09  SYSTEM  PILLING  UNITS  (SPU)  71.36/ 

13.31.09  USER  CPU  TIME  (SFCS)  7.432 

13.31.09  SYSTEM  CPU  TIME  (SECS)  5.733 

13.31.09  USER  MEMORY  USAGE  ( PAGE4SECS  >  1712. C4U 

13.31.09  USER  AVERAGE  WORKING  SET  SIZE  ( PAGES >  22U 

13.31.09  NUMBER  OF  VIRTUAL  SYSTEM  RECUESTS  .90/ 

13.31.09  NUMBER  OF  SMALL  PAGE  FAULTS  850 

13.31.09  NUMBER  OF  DISK  I/O  RECUESTS  27. 

13.31.09  NUMBER  OF  DISK  SECTORS  TRANSFERRED  62v 

13.31.10  *ICOMPLETE»» 
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6. 


I 


Listing  of  Data  File  FINDAT1 


^OLD.FINOAU 
SAMPLE  problem 


16  1 


1 

0 

141 

.76 

1 

2 

IS 

16 

29 

30 

1 

2 

15 

16 

29 

30 

43 

44 

12.00 
0  0 
0.00 
3  4 

17  18 

31  32 

3  4 

17  IS 
31  32 

43  46 


4.71229 

1.00 

20000 

0 

1 

0 

0 

0 

0 

1 

S 

6 

7 

8 

9 

10 

19 

20 

21 

22 

23 

24 

33 

34 

33 

3 

6 

7 

8 

9 

10 

19 

20 

21 

22 

23 

24 

33 

34 

33 

36 

37 

38 

47 

48 

11 

23 

11 

23 

3? 


'V  •. 


12  13 

26  27 

12  13 

26  27 

40  41 


7.  Listing  of  Data  File  FINDAT 


OLDfFINDAT 

LIST 


f 


-8.00 

0.00 

-4.00 

4.00 

0.00 

8.00 

4.00 

0.00 

-8.50 

-11.10 

-32.00 

-32.00 

-20.00 

-20.00 

-8.00 

-8.00 

4.00 

8.00 

8.50 

4.60 

1 

2  4 

4 

7  6 

5 

8  7 

7 

8  10 

7 

12  9 

S 

13  10 

11 

14  16 

12 

13  17 

i  £.  i  ■»  .a 


8.00  -4.00 

-6.00  0.00 

-4.00  4.00 

12.00  11.00 

-12.00 

-32.00  -28.00 

-16.00  -16.00 

-4.00  -4.00 

0.00  4.60 

0.00 

2  3  3 

2  7  4 

3  8  3 

6  9  11 

7  10  12 

11  1 
14  1 

13  1 


4.00 

-8.00 

8.00 

-4.00 

-8.00 

0.00 

8.50 

4.60 

-28.00 

-24.00 

-16.00 

-12.00 

0.00 

0.00 

8.50 

11.10 

1 

4  6 

2 

3  7 

6 

7  9 

9 

12  11 

10 

13  12 

0.00 

8.00 

4.00 

-8.00 

8.00 

-4.00 

0.00 

-4.60 

-24.00 

-24.00 

-12.00 

-8.00 

0.00 

4.00 

12.00 

11.10 

2 


PART  V:  SUMMARY 


56.  A  two-dimensional  finite  element  numerical  simulation  model  (FINITE) 
was  developed  that  calculates  wave  heights  under  combined  refraction  and  dif¬ 
fraction  of  both  long  and  short  waves  approaching  structures  from  any  arbi¬ 
trary  direction.  The  wave  equation  solved  governs  the  propagation  of  periodic, 
small  amplitude  surface  gravity  waves  over  a  variable  depth  seabed  of  mild 
slope.  The  efficient  formulation  of  the  model  permits  the  solution  of  large 
problems  with  relatively  small  time  and  memory  storage  requirements.  A  compu¬ 
tational  scheme  is  employed  that  allows  the  solution  of  practical  problems 
that  typically  require  large  computational  grids. 

57.  Although  the  model  solves  an  equation  that  is  strictly  valid  only 
for  mild  bathymetric  variations,  the  model  can  provide  reasonable  answers  for 
problems  where  there  are  rapid  depth  variations  (at  much  lower  cost  than  re¬ 
quired  by  three-dimensional  models  that  are  appropriate  for  problems  involving 
rapid  depth  variations).  Comparisons  are  presented  between  the  finite  element 
model  calculations  and  an  analytical  solution,  a  two-dimensional  numerical 
solution,  a  three-dimensional  numerical  solution,  and  measurements  from  labora¬ 
tory  studies.  Because  the  finite  element  FINITE  does  not.  provide  a  mechanism 
for  energy  dissipation,  energy  loss  through  wave  breaking  can  be  simulated 
only  by  permitting  waves  to  propagate  out  of  the  computational  region.  The 
program  documentation,  user  guide,  and  sample  problem  output  are  provided. 
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APPENDIX  A:  NOTATION 


a  Wave  amplitude,  ft 

A  Numerical  computational  region,  dimensionless 

A  Wave  amplitude,  it 

Aj  Incident  wave  amplitude,  ft 

B  Numerical  computational  region,  dimensionless 

1/2 

c  Wave  phase  velocity  (g/k  tanh  kh)  ,  ft/sec 

c  Wave  group  velocity  (l/2c(l  +  G))  ,  ft/sec 

8 

F(0)  Mathematical  functional,  dimensionless 

2 

g  Gravitational  constant,  32.174,  ft/sec 

G  Computational  parameter  (2  k  h/sinh  2  kh),  ft/sec 
h  Still-water  depth,  ft 

H  Hankel  function  of  the  first  kind,  dimensionless 
n  1/2 

i  (-1)  ,  dimensionless 

k  Wave  number,  2n/L  ,  1/ft 

k  Incident  wave  number,  2n/L  ,  ft 

o  o 

K  Symmetric  complex  coefficient  matrix,  dimensionless 

L  Arbitrary  wavelength,  ft 

Lq  Incident  wavelength,  ft 

n  Unit  normal,  dimensionless 

nA  Unit  normal  to  the  boundary  separating  computational  regions  A  and 
B  ,  dimensionless 

N  Number  of  node  points  in  the  finite  element  discretization, 
dimensionless 

r  Radial  variables,  ft 

T  Wave  period,  sec 

u  Two-dimensional  velocity  vector,  ft/sec 

or  Constant  coefficients  in  Hankel  functions,  dimensionless 
n 

3  Constant  coefficients  in  Hankel  functions,  dimensionless 
n 

0  Angular  variable  in  polar  coordinates,  deg 

2 

0  Velocity  potential  defined  by  u  =  V0  ,  ft  /sec 

2 

0.  Velocity  potential  in  region  A  ,  ft  /sec 

"  2 

0  Velocity  potential  in  region  B  ,  ft  /sec 

B  2 

0T  Incident  wave  velocity  potential,  ft  /sec 

1  2 
0p  Reflected  wave  velocity  potential,  ft  /sec 


A1 


u)  Angular  frequency,  2n/T  ,  1/sec 

V  Horizontal  gradient  operator,  dimensionless 


A2 
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